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ABSTRACT 


When  the  signal  to  noise  ratio  is  low,  the  detection  of  a  track  in  a  lofargram  by  computer 
is  dinicult.  Actually,  operators  often  face  this  kind  of  problem  when  detecting  an 
acoustic  signature  in  the  noisy  sea  environment.  Furthermore,  operators  must  keep  on 
tracking  lofargram  to  identify  a  given  target.  This  problem  could  be  handled  by  the  au¬ 
tomation  of  lofargram  processing  using  filtering  or  image  processing  techniques.  This 
technique  will  suppress  the  background  and  emphasize  the  spectral  lines  in  the 
lofargram.  Targets  can  be  tracked  by  an  automatic  lofargram  processing  system  up  to  a 
certain  point  at  which  the  system  should  alarm  the  operator.  The  enliancement  proc¬ 
essing  of  lofargrams  using  a  relaxation  method  could  be  one  part  of  the  automatic  sys¬ 
tem  of  lofargram  processing  to  provide  available  target  information  for  good  decisions. 
The  objective  of  this  thesis  is  to  enhance  spectral  lines  of  the  lofargram  by  using  the  re¬ 
laxation  method  which  is  an  iterative  approach  to  line  detection.  This  technique  makes 
probabilistic  decisions  at  every  point  of  the  lofargram  for  each  iteration.  Decisions  are 
adjusted  at  sur'  -'ssive  iterations  based  on  the  results  of  the  previous  iterations.  This 
thesis  tested  algorithms  using  relaxation  method  for  lofargrams.  Some  experimental  re¬ 
sults  of  the  lofargram  processing  are  presented.  Experimental  results  showed  that  an 
edge  relaxation  method  ear.  '.'ield  better  results  than  the  line  relaxation  method.  How¬ 
ever,  a  double  line  detected  from  an  edge  is  still  undesirable  and  requires  future  work. 
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1.  INTRODUCTION 


A.  WHAT  IS  A  LO FAR 

The  term  "LOFAR"  is  flefined  as  "search  technique  using  omnidirectional 
sonobuoys".  It  is  an  acronymn  for  "LOw  Frequency  Analysis  and  Recording."  Every 
ship  including  submarine  create  peculiar  noises  in  the  sea.  The  noise  spectrum  ran  be 
used  to  distinguish  dilTcrent  types  of  snips.  If  we  know  the  noise  frequency  spectrum 
pattern  radiated  from  a  certain  ship,  it  i^  possible  to  identify  this  ship  by  examining  the 
lofargram.  'I'his  is  the  basic  motive  of  using  lofar. 

For  determining  the  characteristic  of  noise  spectra,  a  frequency-time  analyzer  can 
be  used,  which  is  similar  to  those  used  for  speech  analysis.  It  yields  a  plot  of  frequency 
against  time  and  shows  the  intensity  of  the  sound  in  the  bandwidth  of  the  spectrum  by 
darkening  the  record  paper.  This  paper  plot  is  called  a  lofargram. 

Before  the  appearance  of  the  digital  lofargram,  thermal  burning  of  record  paper  was 
used.  But,  the  conventional  paper  gram  was  inconvenient  to  handle  compared  to  digital 
lofargram.  Furthermore,  an  undesirable  smell  was  created  during  operation  because  of 
the  burning  on  the  long  strip  of  paper  using  the  old  method.  Therefore,  the  old  method 
is  almost  replaced  by  the  digital  lofargram.  The  digital  lofargram  is  obtained  by  digitizing 
the  analog  signal  and  displaying  a  gray  .scale  on  the  computer  screen  or  on  a  paper. 

Figure  1  adopted  from  [Ref.  1]  is  an  example  of  a  lofargram  of  the  noise  of  a  large 
surface  ship  obtained  from  a  hydrophone.  As  the  ship  passed  by,  tracks  are  left  in  the 
lofargram.  The  frequency  scale  .-long  the  horizontal  axis  extends  from  0  to  150  Hz  anci 
the  recording  duration  of  the  vertical  axis  was  approximately  1/2  hour.  The  constantly 
spaced  vertical  line  components  marked  by  arrows  aie  blade-rate  lines.  The  linos  marked 
X  are  of  unknown  origin.  Figure  2  adopted  from  [Ref.  2]  indicates  the  block  diagram  of 
a  passive  sonar  to  get  the  lofargram. 

Generally,  the  lofargram  has  characteristics  as  follows 

1.  Signal  to  noise  ratio  is  generally  low. 

2.  A  constant  frequency  tonal  produces  a  darkening  over  time  periods,  which  would 
appear  as  a  vertical  line  on  the  lofargiam  display. 

3.  Multi-lines  could  be  displayed  on  the  lofargram  depending  on  the  spectrum  band¬ 
width. 

4.  Slant  lines  or  slightly  curved  lines  are  due  to  the  Doppler  shift  could  be  shown  on 
the  lofargram. 
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Figure  1.  Lofargrani  of  a  surface  ship  at  a  speed  of  11  knots 


Figure  2.  Passive  sonar  system  block  diagram 


B.  OBJECTIVE  OF  LOFARGRAM  PROCESSING 

When  the  SNR  is  lower  than  3db.  the  interpretation  of  lofargrani  is  difficult.  Actu¬ 
ally,  operators  often  face  this  kind  of  situation  in  the  noisy  sea  environment.  Further¬ 
more,  operators  must  keep  on  tracking  lofargram  to  identify  target.  This  problem  could 
be  handled  by  the  automation  of  lofargram  processing  using  filtering  or  image  process¬ 
ing  techniques  which  can  suppicss  the  background  and  emphasize  the  spectral  lines  in 
the  lofargram.  Target  can  be  tracked  by  an  automatic  lofargram  processing  system  up 


to  a  certain  point  then  the  system  alarms  the  operator.  The  enhancement  processing  of 
a  lofargram  using  the  rela.xation  method  could  be  one  part  of  the  automatic  system  of 
lofargram  processing  to  provide  target  related  information  for  good  decisions. 

The  objective  of  this  thesis  is  to  enhance  the  spectral  lines  of  the  lofargram  by  using 
the  relaxation  method  which  is  an  iterative  approach  to  line  detection.  This  technique 
makes  probabilistic  decisions  at  every  point  of  the  lofargram  for  each  iteration.  Deci¬ 
sions  are  adjusted  at  successive  iterations  based  on  results  of  the  previous  iterations. 

This  thesis  will  present  some  experimental  results  of  the  line  enhanced  lofargrams 
using  the  relaxation  method. 

C.  WHAT  IS  EDGE  AND  LINE  DETECTION 

Edges  are  defined  as  the  discontinuity  of  intensity  in  an  image,  and  arc  basic  parts 
of  the  image  information  in  general.  Both  the  edge  and  line  detection  are  fundamental 
techniques  of  image  processing  which  are  used  to  detect  boundary  of  an  object.  But, 
there  is  an  important  difference  between  edge  and  line  detection.  Edge  detection  aims 
at  finding  the  local  discontinuities  in  an  image.  These  discontinuities  arc  of  interests  be¬ 
cause  they  are  likely  to  occur  at  the  boundaries  of  objects.  The  output  of  an  edge  de¬ 
tection  process  will  be  some  local  line  shown  as  boundaries  of  an  object,  figure  3  shows 
the  step  edge  and  output  profile  of  an  edge  detection  process. 


Figure  3.  Gray  lese!  profile  of  a)  edge  and  b)  line 

In  real  images,  with  noise  and  surface  imperfections,  gaps  between  local  lines  arc  ex¬ 
pected.  Therefore,  a  process  to  organize  the  local  edges  into  aggregates  is  needed  in  a 
line  detection  technique.  Line  detection  is  a  process  where  edges  in  an  image  are  ag¬ 
gregated  to  form  object  boundaries.  These  boundaries  may  not  be  known  a  priori.  But, 
in  many  cases,  they  can  be  appro, ximated  well  by  piecewise  linear  segments.  Mouever, 
it  is  not  feasible  to  simply  fit  linear  segments  to  all  the  edges  in  an  image  and  discard  the 
poor  fits.  It  is  first  necessary  to  aggregate  the  edges  lying  along  a  single  line  or  along 
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another  known  curve.  Proximity  and  the  directions  of  an  edge  or  more  detailed  de¬ 
scriptions  of  the  edges  can  be  used  for  such  aggregation.  Therefore,  edge  detection  and 
line  detection  can  be  put  in  sequence.  Basic  ideas  of  edge  and  line  detection  are  ex¬ 
plained  in  (Ref.  3]. 

Usually  edge  and  line  detection  algorithms  can  be  divided  into  the  following  cate¬ 
gories.  [Ref.  4] 

1.  Detection  of  element. 

Difierential-type  operator. 

Model-fitting  method. 

2.  Enhancement  of  lines. 

Etthancement  in  image  space. 

Enhancement  in  parametric  space  (Hough  transform). 

3.  Connection  of  elements. 

Connections  using  graph  search  techniques. 

Category  1  refers  to  the  process  to  detect  line  or  edge  elements  by  applying  local  oper¬ 
ators  such  as  Roberts  (Ref.  5j  ,  Prewitt  (Ref.  6]  and  Sobel  (Ref.  7J.  An  approximation 
of  the  gradient  for  a  digital  picture  from  Roberts's  opciator  is  given  by 

R{iJ)  =•• 

-  V{«(' +  I)-«(v)}V(g(/,./+  I)-g(/+  lj)f  (1.1) 


In  Figure  4  g{ij)  denotee  the  intensity  of  the  itnage  iof  the  pixel  (i,j),  which  the  direction 
of  the  gradient  is  given  by 


(?  =  _-|-  +  tair'[ 


g{ij+  l)-g(<’+  li/) 
g{i  +lj+  i)-g(v) 


] 


(1.2) 


'fhe  Prewitt  operator  defines  the  magnitude  of  a  gradient  from  Figure  5  by 


Sx  =  (n2+«3-l  fl4)-(«0+«7+<^6) 


(1.4) 


Sy  =  («o+«5+fl4)-(«f0+«l+«7)  (1-5) 

and  the  direction  of  the  edge  is  given  by 
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(1.6) 


.  -I 

6  =  tan  '(  — ) 
->• 


Figure  4.  Robert's  gradient  operator 

The  Sobel  operator  delines  a  gradient  from  figure  5  by 

S  =  I  (c:y+2u,+<3t2)-(Ot,+2as+a4)  | 

+  I  (uo+2a74'a6)-(rt2+2o3-'-c74)  I  (1.7) 

Generally,  a  differential  type  operator  doesn’t  perform  well  compared  to  a  model  fitting 
method.  We  are  going  to  discuss  this  problem  in  the  e.xperiment  part  of  this  thesis. 

Category  2  is  a  smoothing  processes  which  eliminates  noise  components  and  em¬ 
phasize  edges.  In  this  category,  the  iterative  method  performs  edge  enhancement  by  ap¬ 
plying  iterative  processing  based  on  relationships  of  its  neighbour  elements  in  an  image 
space.  Category  3  is  a  process  which  connects  small  line  segments  into  longer  lines. 
Knowledge  or  constraints  related  to  the  sequence  of  elements  to  be  detected  are  used  in 
an  evaluation  function.  Then,  the  connection  of  elements  is  regarded  as  a  graph  search 
problem.  In  this  case,  two  techniques  are  used:  one  directly  proceeds  with  the  line  ex¬ 
traction  while  searching  the  edge  elements  for  brightness  data,  and  the  other  applies  a 
local  edge  operator  to  the  entire  image  before  the  connection  of  elements  is  done. 
[Ref.3] 
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Figure  5.  Apixel  and  Its  eight  neighbors 

In  this  chapter,  the  general  ideas  of  line  and  edge  detection  as  used  in  lofargratn 
processing  arc  explained.  Chapter  II  introduces  the  relaxation  technique,  while  a  modi¬ 
fied  algorithms  using  relaxation  technique  is  introduced  in  Chapter  III.  Experimental 
results  of  lofargram  processing  are  given  in  Chapter  IV.  Chapter  V  is  the  conclusion 

of  this  thesis. 
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II.  RELAXATION  METHOD 


A.  GENERAL  BACKGROUND 

Relaxation  tccliniqucs  were  first  used  by  Waltz  [Ref.  8]  for  the  description  of  solids 
and  were  expanded  by  Zuckcr  [Ref.  9]  for  a  variety  of  applications.  Relaxation  is  an  it¬ 
erative  approach  to  assign  each  pixel  to  categories  by  assigning  the  neighboring  pixels 
in  a  "compatible"  way.  Ihis  means  that  each  pixel  has  interaction  with  its  neighboring 
pixels,  and  the  degree  of  interaction  is  considered  in  terms  of  the  compatibility  of  each 
pixel  to  its  neighboring  pixels,  fhis  method  takes  a  set  of  probabilities  for  each  pixel  to 
belong  to  a  possible  class,  and  uses  an  iteiative  technique  to  update  the  probabilities. 

An  important  clement  of  the  method  is  a  set  of  compatibility  measures  c(/<,.s;/,/)  that 
gives  the  compatibility  of  assigning  pixel  k  to  class  s  and  assigning  the  neighboring  pixel 
/  to  class  /.  It  is  assumed  that  c(/<,s;/,/)  is  in  the  range  (-1,+  1)  where  -1  implies  a  strong 
incompatibility,  +  1  implies  a  strong  compatibility,  and  0  implies  neutrality.  p'^{k,s)  is  an 
initial  estimate  of  the  probability  that  pixel  k  belongs  to  class  s.  Let  q'{k,3)  is  the  incre¬ 
ment  of  the  probability  and  N  bt  the  number  of  neighbors  for  each  pixel;  a  neighbor¬ 
hood  of  3x3  or  5x5  is  usually  employed.  Then,  a  sequence  of  estimates  p'{k,s)  is 
computed  iteratively  as  follows 

1.  For  each  pixel  /<,  compute  the  "neighbor  compatibility"  between  pixel  k  t.nd  its 
neighbor  pixel  /  in  class  s,  which  results  in  the  probability  of  pixel  k  in  class  s. 

=  (2-1) 
/  f 

where  the  summation  is  over  all  possible  neighbor  pixels  /  and  all  classes  t.  If  only 
two  classes(for  example,  edge  and  non-edge)  are  considered,  then  there  are  only 
four  cases  of  neighbor  compatibility  c{k,e\l,e)  means  the  compatibility  of  an  edge 
pixel  k  being  a  neighbor  of  an  edge  pixel  /.  c{k,c:,l,n)  means  the  compatibility  of  an 
edge  pixel  k  being  a  neighbor  of  a  non-edge  pixel  I .  c{k,n;l,e)  means  the  compat¬ 
ibility  of  a  non-edge  pixel  k  being  a  neighbor  of  an  edge  pixel  / .  c{k,iv,l,n)  means 
the  compatibility  of  a  non-edge  pixel  /<  being  a  neighbor  of  a  non-edge  pixel  /  . 
Therfore,  equation  (2.1)  can  be  expanded  as  follows 
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(2.2) 


q\k,e)  =  —  ^c(/c, <?;/,«)//(/, e)+c(/<,e;/,H)/?V.«) 
/ 


q'{k,n)  =  ~  '^c{k,n;l,e)p\l,e)+c{k,n\l,n)p{l,it)  (2.3) 

l 


2.  Update  the  probabilities  p'*'{k,s) 

r+\,,  _  p\k,s)ilWil<,s)) 

P  \ktS) 

2_^p\k,s)[l+q\k,s)) 

S 


For  two  classes; 

r+1.^  .  ^ _ />V<,g)[l+?Vc,c)] _ 

^  />V<,e)L  •+'/'(/<.«)]+/> '(^.«)[J+<7'(^.»)] 


(2.4) 


(2.5) 


Equauon  (2.5)  gives  the  probability  that  pixel  k  is  an  edge.  Then  the  probability 
that  pixel  /<  is  a  non-edge  is  defined  by 

//+'(/t,/0=l-/^'(A,e)  (2.6) 


As  the  process  is  iterated,  these  probabilities  must  stay  in  the  range  from  0  to  1.  The 
formulas  given  above  are  the  basic  expression  of  the  relaxation  technique.  This  concept 
can  be  expanded  and  modified  for  other  applications.  This  technique  has  the  following 
characteristics: 

1.  High  compatibility  pixels  tend  to  reinforce  each  other  and  the  low  compatibility 
pixels  tend  to  discourage  each  other. 

2.  The  degree  of  reinforcing  oi  discouraging  done  by  a  neighbor  is  proportional  to  its 
own  probability  of  assignment  of  that  class. 

3.  The  probabilities  p{k,s)  remain  in  the  range  from  o  to  1  and  sum  to  1  over  all  of 
class  s. 

B.  APPLICATIONS  OF  RELAXATION  METHOD 
I.  Line  enhancement 

One  of  the  applications  in  which  relaxation  techniques  liave  been  proven  useful 
is  the  detection  of  long  smooth  curves  in  an  image.  Initially,  line  detection  operators 
are  applied  to  the  image,  and  their  outputs  are  used  to  determine  an  initial  probability 
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of  the  point  lying  on  a  curve  with  a  given  orientation.  Tliese  probabilities  are  then  iter¬ 
atively  reinforced:  the  probability  of  a  point  lying  on  a  curve  is  reinforced  by  the  prob¬ 
abilities  of  the  other  points  lying  on  curves  that  smoothly  continue  it.  After  iterations 
of  the  procedure,  points  that  lie  on  smooth  curves  tend  to  have  high  probabilities  of 
being  curve  points,  while  the  other  points  do  not. 

An  application  of  the  relaxation  technique  to  the  detection  of  major  edges  in 
an  image  will  be  presented  in  the  next  chapter  in  detail.  Initially,  a  gradient  edge  opera¬ 
tor  is  applied  to  the  image.  This  provides  information  about  edge  strength  and  orien¬ 
tation  at  each  point.  These  probabilities  arc  then  iteratively  reinforced. 

The  general  idea  of  the  edge  reinforcement  process  is  as  follows. 

a.  Magnitude  and  direction 

First,  the  magnitude  and  direction  of  the  image  gradient  are  computed  for 
each  pixel  /<.  The  magnitude  of  pixel  /c,  divided  by  the  maximum  of  the  magnitudes  over 
the  entire  image,  defines  the  initial  probability  p^{k,s)  of  a  pixel  k  being  an  edge. 

b.  Reinforcement  of  edges 

The  reinforcement  process  defines  a  new  edge  probability  of  a  pixel  k  in 
terms  of  the  old  probabilities  at  the  pixel  k  and  its  ncignbor  pixels.  The  computation  of 
the  new  edge  probability  can  be  broken  into  steps  as  follows. 

(])  Interactions  between  edge  and  non-edge.  The  interaction  between 
edge  and  edge,  edge  and  non-edge,  non-edge  and  edge,  non-edge  and  non-edge  depend 
on  the  edge  direction  at  the  two  positions,  on  the  orientations  of  the  line  joining  the 
points,  and  on  the  distance  between  them. 

(2)  Computation  of  the  new  probability.  To  compute  the  new  edge  prob¬ 
ability  obtained  from  the  iterations,  first,  compute  weighted  sums  of  the  edge  and  non- 
edge  probability  increments,  q'{k,e)  and  q'{k,n).  These  sums  .arc  then  normalized  to  lie 
in  the  rangc(-l,+  1)  according  to  (2.4).  They  are  used  to  update  the  edge  and  non-edge 
probabilities,  while  the  updated  values  arc  normalized  so  that  they  can  sum  to  1. 

The  real  issue  in  this  thesis  is  to  search  for  good  line  detection  tech¬ 
niques  applicable  to  lolargrams,and  not  to  edge  detection.  Therefore,  we  are  going  to 
concentrate  on  the  reinforcement  of  the  relaxation  method  as  applied  to  lofargratn. 

2.  Labeling 

There  is  another  application  using  the  relaxation  technique,  which  is  called  "la¬ 
beling."  The  objective  is  to  assign  labels  to  objects  in  an  image.  [Ref.  10] 
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a.  Components  of  labeling 

In  a  finite  set  of  relations  between  objects,  the  objects  usually  correspond 
to  entities  to  be  labeled.  The  objects  are  often  geometrically  or  topologically  related  to 
each  other.  An  input  scene  is  thus  a  relational  structure  of  all  objects.  In  the  simplest 
case,  each  object  is  to  be  assigned  with  a  single  label.  Labels  may  be  weighted  with 
"probabilities"  indicating  the  "probability  of  an  object  having  that  label.  "Constraints 
determine  what  labels  may  be  assigned  to  an  object.  A  basic  labeling  problem  is  then: 
Given  a  finite  input  scene  (relational  structure  of  objects),  a  set  of  labels,  and  a  set  of 
constraints,  find  a  consistent  labeling. 

There  are  several  types  of  labeling  using  relaxation. 

b.  Discrete  labeling 

It  is  a  parallel  iterative  algorithm  adjusting  all  object  labels.  All  possible 
labels  to  each  object  arc  assigned  in  accordance  with  constraints.  Iterations  are  per¬ 
formed  until  a  globally  consistent  labeling  is  found.  In  parallel,  from  each  object's  label 
set  all  the  labels  that  are  inconsistent  with  the  current  labels  of  the  rest  of  the  relational 
stiuctuie  are  climinatecl. 

c.  Linear  labeling 

The  labeling  process  starts  with  an  initial  assignment  of  weights  to  all  ob¬ 
jects.  Weight  are  reminiscent  of  probabilities,  reflecting  the  "probability  that  a  label  is 
correct."  ilere  p  refers  to  probabilily-like  weight  rather  than  to  the  value  of  a  probability 
density  function.  Let  a  relational  structure  with  n  objects  be  given  by  <?* ,  /c  =  each 
with  nt  discrete  labels  and  let  p^{?.)  denote  the  weight,  or  the  "probability"  that  the 
label  is  correct  for  the  object  a*  and  restricted  as  follows 

0  ^  I  (2.10) 

(211) 

.1 

The  linear  labeling  operator  is  based  on  the  compatibilities  of  the  labels, 
which  serve  as  the  constraints,  A  compatibility  p^|  looks  like  a  conditional  probability. 

(2-12) 

A 
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The  pyf^k  I  >T)  may  be  interpreted  as  the  conditional  probability  that  object  o*  has  label 
X  given  that  another  object  a,  has  label  X\  The  operator  iteratively  adjusts  label  weights 
in  accordance  with  other  weights  and  the  compatibilities.  A  new  weight  p^iX)  is  com¬ 
puted  from  old  weights  and  compatibilities  as  follows. 

(2-13) 

/ 

The  c„  are  coefficient  such  that 


d.  Nonlinear  labeling 

Ii  the  compatibilities  are  allowed  to  take  on  both  positive  and  negative 
values,  we  can  express  strong  incompatibility  and  strong  compatibility.  Denote  the 
compatibility  of  the  event  "label  X  on  aC  with  event  "label  X'  on  a"  by  r^|{X,  >!')•  the 
two  events  occur  together  ol'ten,  r*,  should  be  positive.  If  they  occur  together  rarely, 
should  be  negative.  If  they  are  independent,  r*,  should  be  0.  ’I'he  compatibilities  are  based 
on  correlations 


(2.15) 

a'(A')  =  /;(A-,A-)-(p(A0)' 

(2.16) 

(2.17) 

The  r^|  is  used  to  obtain  the  positive  or  negative  change  in  weight. 

/  A- 

The  q\{X)  is  the  increment  of  the  weight.  The  probability  change  is  as  follows 
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(2.19) 


Basically,  labeling  is  dinbrcnt  with  the  detection  problem.  Labeling  is  a 
process  to  assign  the  correct  label  to  the  correct  position  while  detection  is  a  process  to 
take  the  desired  portion  of  a  scene  having  a  mixed  feature.  Therefore,  labeling  is  a  dif¬ 
ferent  technique  for  line  detection. 

The  general  concept  of  the  relaxation  technique  and  its  applications  were 
discussed  in  this  chapter.  In  the  next  chapter,  line  and  edge  detection  algorithms  using 
relaxation  technique  will  be  discussed  in  details. 
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III.  RELAXATION  ALGORITHMS  APPLIED  TO  LINE  DETECTION 


There  are  many  algorithms  already  developed  for  line  detection,  in  this  thesis,  the 
objective  ij  to  select  and  test  some  of  them  for  the  lofargram  processing. 

A.  EDGE  ENHANCEMENT(RXEG) 

This  algorithm  is  derived  from  the  original  idea  developed  by  Schachter,  et  al  to  re¬ 
inforce  the  contin.mus  edges  detected  by  a  gradient  operation.  RXEG  is  an  acronymn 
for  "Relaxation  method  for  EdGe  dctcction".Edgcs  reinforce  other  edges  that  interact 
with  nearby  non-edge  points  in  specified  ways.  The  gradient  of  the  edges  are  also  iter¬ 
atively  adjusted.  (Ref.  11) 

I.  Initial  edge  probability 

A  digital  gradient  operation  is  applied  to  the  given  image  f.  if  we  denote  the  .t 
and  y  components  of  the  gradient  by  A/ and  A/,  then  the  magnitude  and  direction  of 
the  gradient  .ire  given  by 

(3.1) 


(3.2) 

'1  he  mag  of  die  gradient  indicates  the  strength  of  the  edge.  1  he  angle  0  of  the  gradient 
is  perpendicular  to  the  edge  direction.  Edge  direction  is  obtained  by  adding  90'  to  the  0 
angle  of  the  gradient. 

We  defined  the  "probability"  of  an  edge  at  a  given  point  k  by 


p{k,e) 


mag{k) 

iriii\{mag{l)) 


(3.3) 


where  the  max  is  taken  over  the  entire  image.  The  probability  of  a  nonedge  at  a  pixel  k 
is  defined  as  p{k,n)  —  \—p[k,e). 

2.  Pixel  and  neighbor  interaclion 

The  centered  pixel  k  could  have  four  kinds  of  interactions  with  its  neighbor  pixel 
/  such  as  edge/edge,  edge/non-edge,  ncn-cdgc/cdge  and  non-edge/iion-edgc. 
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0.  Edgejcdge  inteiaction 

k  is  the  centered  pixel  position,  and  /  is  one  of  the  neighbor  pixel  position. 
Let  a  be  the  edge  direction  at  k,  fi  the  edge  direction  at  /,  y  be  the  direction  of  the  line 
joining  A  to  /,  D  the  chessboard  distance  from  k  to  /,  i.e,  max(  1  k^-l,  | ,  |  l<y~ly  1 ).  Then 
the  edge/edge  reinforcement  process  between  the  points  k  and  /  has  strength  given  by 

c{k,e\l,e)  =  cos(a-y)  cos(/?-y)/2”  (3.4) 

To  sec  the  significance  of  this  definition,  a  few  simple  examples  are  consid¬ 
ered.  In  these  examples,  the  arrows  indicate  the  direction  along  the  edge,  with  the  dark 
side  of  the  edge  on  the  left.  These  examples  showed  that  parallel  and  perpendicular 
edges  have  no  effect  on  one  another. 


Case 

a 

P 

y 

cos(a  -  i)  cos(p  -  Y  )/2‘^ 

n 

90 

90 

0 

0 

90 

270 

0 

0 

A 

A 

90 

90 

90 

1/2^' 

A 

90 

270 

90 

-1/  2^ 

90 

0 

0 

0 

■■ 

90 

180 

0 

0 

Table  1.  Examples  of  edge  /  edge  interaction  [Ref.  8] 


b.  Edgej non-edge  interaction 

Besides  the  edge/edge  interaction  which  occurs  between  the  edge  probabili¬ 
ties  p{k,e)  and  j){l,e)  there  are  also  interactions  involving  the  non-edge  probabilities 
p{l,n) .  The  edge  probability  at  pixel  k  is  weakened  by  the  non-edge  probability  at  pixel 
/  to  the  degree  c{k,e\l,n)  defined  by 

c{k,e-,l,n)  =  min[0,  -  cos(2a-2y)/2^]  (3.5) 
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In  the  examples,  non-edge  neighboring  points  collincar  with  the  center  edge  points 
weaken  them,  whereas  non-edge  points  alongside  edge  points  have  no  eflect  on  them. 


Case 

a 

90 

270 

• 

90 

A 

90 

-  cos(2a  “  2y)  /2^ 

c(k,  e;  1,  n) 

H2^ 

0 

112^ 

0 

-\I2^ 

-1/2*^ 

-1/2^ 

-112^ 

Tabic  2.  Examples  of  edge  /  non-edge  interaction  (Ref.  8] 
c.  Non-edgejedge  hucraction 

'1  he  non-edge  probability  at  pixel  k  is  alTectcd  by  the  edge  probability  at  the 
neighbor  pixel  /  to  the  degree  c[k,n\l,e)  defined  by 

./  .  (l-cos(2/?-2r)) 


Case 


■A 


A 


( 1  -cos(2a  -  2y))  /2 


Tabic  3.  Examples  of  non-edge  /  edge  interaction  (Ref.  8] 


I'rom  the  examples,  neighbor  edge  points  alongside  the  center  point  strengthen  them, 
While  edge  points  collincar  with  non-edge  points  have  no  elfect  on  them. 
d.  Non-edge  I  non-edge  interaction 

1  he  non-edge  probabilities  at  centered  pixel  k  and  the  neighbor  pixel  I  re¬ 
inforce  each  other  to  the  degree  c{k,n',l,n)  defined  by 


c{k,ir,l,n)  =  ~ 
2 
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Non-edge  probabilities  arc  reinforced  by  the  otlier  nearby  non-edge  probabilities. 
c(/-, «;/,«)  is  directionally  independent. 

3.  Combined  reinforcement  proce.ss 

For  each  pixel  k,  the  net  elfcct  of  its  neighboring  pi,xcl  on  its  edge  probability 
p{k,e)  and  non-edge  probability  pik,n)  =  is  computed  as  follows 


q'{k,e)  = 

/ 


+y  ^C2p'il,it)c{k,c;l,it) 


(3.8) 


=  Y^Qp\l,<!)c{k,ivJ,e) 

I 


I 


where  C,,  C\,  Cj,  C*  are  constants  whose  sum  is  equal  one.  The  standaid  values  of 
C,  =  0.S6()  ,  C2  =  0.I24,  0,  =  0.005,  and  C«  =  0.005  was  used  in  the  paper  by  Bruce 
J.Schacliter,  et  al  [Ref.  11].  The  results  of  the  iteration  process  are  somewhat  sensitive 
to  the  choice  of  the  Cs.  For  example,  if  C,  is  too  large,  the  edge  will  thicken  and  will 
be  extended  into  non-edge  points;  while  if  C4  is  too  large,  gaps  will  appear  at  W'eak  spots 
in  the  edges  and  at  sharp  angles. 


q\k,e)=-r  - 

^  I  q\k,c)  I  -t- 1  q\k,n)  I 

I  q\Ke)  I  -1- 1  q\k,n)  \ 
p'{k,e)  =p\k,e)l\-\-q'{k,e)'] 
p'{k,n)==  p\k,ii)iU-q’{k,n)] 


(3.10) 

(3.11) 

(3.12) 

(3.13) 
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(3.14) 


p'[k,e)+p'(k,it) 


This  process  is  then  iterated  by  using  p'^'{k,e)  iti  place  of  p'{k,e),  and  by  using 
in  place  of  p'{k,n).  Wc  also  compute  the  estimated  edge  direction  A',(A)  and 
at  each  point 

A'#)=liy(/.g)cos(0(/<)) 


+^p\l,e)c\k,c\l,e)  cos(0(/)) 


/ 


A'^,(/c)  =  ]Vp\l,e)  %m[0{k)) 


(3.15) 


+^//(/.c)c(/v,<?;/,r)  sin(0(0)  (3.16) 

I 

r'(A)  =  tan-'(A',(.v^')/A>j-))  (3.17) 

For  large  values  of  the  constant  W,  0'^'  is  close  to  0;  while  for  small  values,  it  is  strongly 
influenced  by  the  neighboring  O's. 

n.  LINE  AND  CURVE  ENIIANCEMENTlRXLN) 

In  section  A,  only  two  categories  arc  associated  with  the  image  pixel,  (i.e,  edge  and 
no.i-cdge).  In  this  section,  more  categories  arc  considered  in  the  RXLN  relaxation  pro¬ 
cedures.  RXLN  is  an  acronymn  for  "Relaxation  method  for  LiNc  detection." 

This  algorithm  was  devlopcd  by  Zuckcr,  Ilununcl  and  Roscnfeld  (Ref.  I2J.  The  re- 
la.xation  process  i>  applied  to  the  detection  of  .smooth  lines  and  curves  in  noisy  images. 
Nine  class  labels  associated  with  each  image  point  will  be  considered.  Light  classes  in¬ 
dicate  lines  at  various  orientations  and  one  indicates  the  no-Iine  case.  In  the  relaxation 
process,  interaction  takes  place  between  the  probabilities  at  neighboring  points.  This 
permits  line  segment  with  compatible  orientations  to  strengthen  one  another.  Similarily, 
no-line  class  is  reinforced  by  neighboring  no-line  class  and  weakened  by  the  neighboring 
oriented  line  classes. 
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1.  Initial  piobabilities 

At  each  point,  probabilities  are  assigned  to  the  nine  classes  which  are  lines  in 
eight  possible  directions  and  a  no-linc  class.  The  initial  probability  for  each  class  is  ob¬ 
tained  by  evaluatng  a  nonlinear  line  detector  as  shown  in  Figure  6  at  every  picture  point 
along  the  eight  orientations.  Eight  sets  of  pixels  representing  eight  orientations  are  se¬ 
lected  and  shown  in  Figure  6.  Each  set  of  pixels  arc  rearranged  as  the  template  shown 
in  Figure  6  (b).  When  the  condition  in  Figure  6  do  not  hold,  the  response  is  zero.  When 
they  do  hold,  the  response  is  calculated  by 

R=={D+E+H)-~{A+D+G+C+F+r)  (3.18) 

At  each  pixel,  the  detector's  response  are  computed  for  every  orientation,  'i  hen  only  one 
orientation  which  has  the  maximum  detector's  response  is  selected  as  the  orientation  of 


the  pixel.  [Ref.  13] 

.Temalala 

Conditions 

A 

B  C 

A 

<  B  >  C 

D 

E  F 

D 

<  E  >  F 

G 

H  1 

(a) 

G 

<  H  >  1 

B  B 

B 

B 

H 

H 

E  E 

E 

HE  H 

E  B 

HE  E 

E 

H  H 

H 

B  B 

B 

(1)  (2) 

(3) 

(4) 

(b) 

(5) 

(6)  (7) 

(8) 

Figure  6.  Nonlinear  line  detector,  (a)  Nonlinear  line  dctecioi  for  the  vertical  ori¬ 
entation.  When  the  conditions  do  not  hold,  the  response  is  zero;  when 
they  do  hold,  the  response  is  calculated  by  (3.18).  (b)  Eight  orientations 
of  the  detector. 
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To  describe  this  process,  let  d,(k)  denote  the  output  of  the  nonlinear  detector  for 
orientations  s  at  pixel  k  .  Let  p{k,s)  denote  the  probability  that  a  pixel  k  belongs  to  class 
s.  Then  the  initial  probabilities  can  be  obtained  by  scaling  the  detector's  responses  at 
each  position  by  the  maximum  over  the  entire  picture.  The  probability  that  pixel  k  has 
class  s  is  expressed  as  follows 


p{k,s) 


nvxMk)) 

j=i _ 

max((max)(/j(A)) 

S=\  It 


8 

1=1 


(3.19) 


max(rf,(-^))  means  the  one  selected  maximum  value  over  all  orientations.  max(<y,(/0) 

j-i  * 

means  the  maximum  value  of  oriented  s  over  all  image  pixels,  and  the  probability  that 
pixel  k  has  no-line  is  defined  by 


8 


p{k,9)=  ]-2_f{k,s) 


1=1 


(3.20) 


2.  Conipalibility  coefficients 

To  update  the  probabilities,  the  compatibility  relation  between  a  pixel  and  its 
neighbors  must  be  specified..  The  compatibilities  depend  only  on  the  relative  orientations 
of  the  neighboring  point.  These  compatibilities  can  be  specified  in  the  following  way.  If 
two  neighboring  line  segments  arc  oriented  in  the  same  direction  or  close  to  the  same 
direction,  they  add  support  to  one  another.  If,  on  the  other  hand,  two  neighboring  seg¬ 
ments  arc  oriented  perpendicularly  to  one  another,  they  will  reduce  support  to  each 
other.  All  other  pairs  of  line  segments  are  distributed  between  these  two  extremes.  In  this 
algorithm,  there  arc  three  types  of  compatibility  coefficients. 
a.  Compatibility  coefficients  between  lines 

Figure  7  shows  the  compatibility  coefficients  between  the  lines  in  five  ge¬ 
ometrical  relationships.  As  standard  values,  from  left  to  right:  1.0,  0.5,  0.05,  -0.15,  -0.25 
are  used.  When  straight  line  enhancement  is  desired,  1.0,  0.0,  -0.1,  -0.17,  -0.22  ate  used. 
[Ref.  13] 
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b.  Compatibility  coefficients  between  fine  and  no-lines 

As  its  standard  value,  -0.1  is  used.  Decreasing  the  value  to  -0.5  increases  the 
elTect  to  fill  the  gap  between  line  segments.  [Ref.  13]  The  standard  value  is  0.25. 

I  /  / 

'■0  .5  ,05  -.15  -.25 

Figure  7.  Compatibility  weights  between  line  labels 

f.  Compatibility  coefficients  between  no-lines 

Increasing  the  value  up  to  0.5  improve  the  noise  efiect  but,  may  cause  era¬ 
sure  of  line  labels.  [Ref  12] 

For  each  pixel  /<,  its  neighbor  pixel  has  9x9  cases  of  interaction  with  the 
center  pixel.  If  neighborhood  size  is  5.x5.  ’Ihcrc  are  5x5x9x9  cases  of  compatibility  coef¬ 
ficients, 

c{k,\\\,\)  .  .  c(A-,l;25,l) 

t  •  •  • 

•  •  *  • 

c(/<,9;l,9)  .  .  c(*,9;25,9) 

3.  Updating  probability 

The  updating  process  can  be  expressed  in  terms  of  compatibility  functions.  Let 
q'{k,s)  represents  the  increment  applied  to  p'{k,s),  which  is  the  probability  that  pixel  k 
belongs  to  class  s.  The  probabilities  of  each  pixel  in  each  class  are  computed  as  follows 

5;/, /)//(/,/)  (3.21) 

/  t 


where  s  repr^'sents  nine  classes  of  reference  pixel  and  t  represents  nine  classes  of  neigh¬ 
boring  pixels  and 


(3.22) 


r+V/  s  p\k,s)li-\-q\k,s)'] 

^p\k,s)[l+q\k,s)^ 

S 

C.  LINE  AND  CURVE  ENHANCEMENT(RXLA1) 

This  algorithm  is  delivered  from  the  original  idea  developed  by  S.Peleg  and 
A.Rosenfeld  (Ref.  14).  RXLAl  is  an  acronynin  for  "Relaxation  method  for  Line  And 
curve."  Basic  concept  using  in  this  algorithm  is  the  same  as  that  used  in  the  RXLN(sce 
section  B)  except  that  the  method  for  compatibility  cociricient  is  different. 

1.  Compatibility  coefficient 

One  possible  interpretation  of  the  compatibilies  is  in  terms  of  statistical  corre¬ 
lation.  Couclation  has  properties  as  follows 

a.  If  class  s  and  class  t  arc  compatible  for  pixel  k  and  pixel  /,  c{k,s\l,t)  >  0 

b.  If  class  s  and  class  t  arc  incompatible  for  pixel  k  and  pixel  /,  c{k,s\l,t)  <  0 

c.  If  neither  class  is  constrained  by  the  other,  c{k.s',t,t)  =  0 

d.  The  magnitude  of  c{k,s',l,i)  reptesents  the  strength  of  the  compatibility. 

Estimates  of  the  correlation  cocllicicnts  derived  from  analyzing  the  initial 
probability  are 


c{k,s\l,t) 


[/>(/c,.v)-/T(/<,.s)3[/>(/,0-p(/.0] 

(7(S)0(/) 


(3.23) 


/'(/{, s)  is  the  initial  estimate  of  the  probability  of  pixel  k  with  class  s,  p{k,s)  is  the  average 
of  p(/<,s)  for  all  pixels  k,  a(.s)  is  the  standard  deviation  of  p{k,s).  If  the  neighborhood  size 
is  5  X  5,  5  X  5  X  9  X  9  cases  exist  for  each  pixel  as  below 

c(/<,I;l,l)  .  .  c(/v,l;25,l) 

•  »  •  • 

•  •  •  • 

c(^,9;l,9)  .  .  c(^,9;25,9) 

Then,  for  every  case  and  for  all  pixel  /<  we  have  the  summation  as  given  by 
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c(A,l;I,l)  .  .  c(/<,l;25,l) 

•  •  •  • 

c(A„l;l,l)  .  .  c(/<„l,25,l) 

•  •  *  * 

•  •  •  • 

C(A,9;1,9)  .  .  c(/<,9;25,9) 

•  •  •  « 

c(A„9;1,9)  .  .  c(A„9;25,9) 

c(A-2,1;1,1)  .  .  c(A2,1;25,1) 


+ 


c(/<2.9;1,9)  .  .  c(/<2.9;25,9) 


+  •  •  • 


Thcrerorc 


k  _  _ 


(7(5)a(f) 


(3.24) 


In  tins  algorithm,  there  is  the  cllect  of  dominance,  that  is,  wlicn  one  class 
dominates  the  picture,  its  correlation  coellicicnts  with  all  other  classes  arc  high.  Thus  in 
the  relaxation  updating,  the  dominant  class  gets  most  of  the  support,  and  after  a  few 
iterations  almost  all  the  pixels  have  this  class.  '1  he  clfect  of  dominance  among  classes 
can  be  alleviated  by  weighting  the  c{k,s]l,i)  by  the  probabilities  that  the  corresponding 
classes  do  not  occur.  This  greatly  reduces  the  values  of  the  coellicicnts  involving  domi¬ 
nant  classes,  but  will  only  slightly  reduce  the  values  of  the  coellicicnts  involving  rare 
classes.  '1  he  compatibility  coellicicnts  used  in  this  algoritm  arc 

c(A,s;/,0  =  C1-p(/<,5)][1-/T(/,/)]  X  c{k,l,s,t)  (3.25) 


Two  line  and  edge  detection  algorithms  using  relaxation  technique  were  dis¬ 
cussed  in  this  chapter.  These  algorithms  will  be  applied  to  the  lofargram  for  line  de¬ 
tection.  Experimental  results  of  these  algorithms  will  be  discussed  in  Chapter  IV. 
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IV.  EXPERIMENTAL  RESULTS  AND  A  MODIFIED  RELAXATION 

APPROACH 


In  the  experiments,  four  difFerent  relaxation  algorithms  are  applied  to  lofargrams 
having  difFerent  S/N.  Both  real  and  artificially  created  lofargrams  were  used  as  test  im¬ 
ages.  As  discussed  previously,  two  algorithms  RXLN  and  RXLAl  arc  line  detection  al¬ 
gorithms  while  the  third  algorithm  RX12G  is  an  edge  detection  algorithm.  1  he  purpose 
of  the  experiment  is  to  test  these  three  algorithms  to  sec  how  they  work  for  lofargrams. 
■fhe  secondary  purpose  of  these  experiments  is  to  obtain  thin  lines  for  loFargram  inter¬ 
pretation.  Because  the  results  of  these  algorithms  were  not  satisfactory',  improvement 
oFthc  most  suitable  algorithms  arc  suggested. 

The  results  of  the  line  detection  algotithms  RXLN  and  RXLAl  were  not  encour¬ 
aging  when  the  S,'N  is  low.  'Fhe  reason  will  be  explained  later  in  this  chapter.  Fhe  results 
improve  as  the  S/N  increases.  One  effect  is  that  the  noise  was  removed  very  slowly  in 
iterations  of  reinforcement.  On  the  other  hand,  the  edge  detection  algorithm,  RXHG 
showed  better  results  than  the  above  two  algorithms  in  the  presence  of  low  S/N.  But  the 
results  were  still  not  ideal.  A  double  line  was  shown  when  the  S/N  was  low’  because 
RXLG  was  the  edge  detection  algorithm.  'Fheieforc,  some  modification  is  needed  for  line 
detection  in  tlie  RXIiG.  This  will  be  di.scusscd  in  details  later. 

Two  parts  in  the  RXIIG  algorithm  were  modified.  The  first  one  is  the  initial  magni¬ 
tude  of  the  pixel,  since  the  outputs  arc  greatly  influenced  by  the  initial  magnitude  of  the 
input  image.  I  hc  other  is  the  initial  direction  of  the  pixel.  Because  RXIIG  uses  a  difl'cr- 
ential  type  of  edge  detection  operator,  the  initial  direction  in  RXEG  is  not  appropriate 
for  the  lofargrarn  line  detection.  In  RXIIG  the  edge  direction  has  a  dark  side  which  in¬ 
dicates  the  outside  of  the  object  and  a  blight  side  which  indicates  the  inside  of  the  object 
body.  But  the  direction  of  a  line  simply  indicates  the  potential  direction  of  a  track  in  a 
lofargrarn.  Detail  will  be  discussed  in  section  C. 

In  a  search  for  a  better  initial  magnitude  and  direction,  many  other  algorithms  w'ere 
tested.  It  assures  that  the  magnitude  and  diicction  output  from  EGPR  which  is  another 
line  detection  can  be  useful.  The  RXEG  was  modified  and  tested  accordingly. 

A.  LINE  DETECTION  ALGORITHMS  (RXLN,  RXLAl) 

Both  RXLN  and  RXLAl  are  line  detection  algorithms  which  use  similar  methods 
for  initial  probability  and  compatibility  coefficient  calculations  as  explained  in  Chapter 
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HI.  But  the  applications  of  these  are  very  discouraging.  Figure  8  shows  a  real  lofargram 
as  the  input  image  and  its  initial  magnitude  calculated  in  RXLN.  Figure  9  shows  the 
results  of  two  iterations  of  reinforcement.  And  Figure  10  shows  an  artificial  lofargram 
with  S/N  of  6  db  and  its  initial  magnitude  image.  Figure  1 1  show  the  results  of  the  four 
iterations  of  reinforcement.  And  Figure  12  show  the  results  for  3  db  and  6  db  of  S/N. 

The  initial  magnitude  of  input  image  shows  that  the  result  is  vciy  noisy.  In  these 
algoritlims,  there  is  the  eficct  of  dominance.  When  one  class  dominates  the  picture,  high 
value  is  assigned  to  the  dominating  class  at  each  pixel  and  low  values  are  assigned  to  the 
other  classes.  Thus,  the  dominant  class  gets  the  most  support  in  the  relaxation  updating 
process.  And  after  a  few  iterations  almost  all  the  pixels  have  the  same  class.  Figure  8 
shows  tliat  the  noise  dominates  the  picture  when  the  lofargram  has  low  S/N.  Figure  10 
(b)  still  shows  that  the  noise  dominates  the  picture.  When  the  S/N  is  higher,  the  signal 
showed  up  weakly  as  shown  in  Figure  1 1.  But  noise  still  dominates  the  picture.  Gener¬ 
ally  noise  dominates  the  lofargram.  Noise  sttongly  dominates  the  lofargram  when  the 
S/N  is  low.  Hence  the  RXLN  and  RXLAl  algorithms  are  not  suitable  for  lofargram 
processing. 

B.  THE  RESULTS  OF  THE  EDGE  DETECTION  ALGORITHM  RXEG 

Figure  13  shows  the  initial  magnitude  of  the  RXBG  of  the  input  image.  Figure  14 
shows  the  results  of  the  reinforcement  applied  to  Figure  13.  In  this  implementation, 
neighborhood  size  is  5x5,  then  the  24  neighboring  pixels  inlluenced  the  centered  pixel. 
Figure  15  shows  the  artificial  lofargram  input  image  and  its  initial  magnitude  with  6db 
of  S/N.  The  iteration  scheme  has  considerable  noise  removal  power  as  shown  in  the 
example  of  Figure  14.  Figure  16  show's  that  the  results  of  the  artificial  lofargram  with 
higher  S/N  can  yield  a  stronger  response.  A  double  line  in  Figure  14  (b)  and  Figure  16 
indicate  that  the  edge  w'as  detected  from  the  input  image.  Because  RXFG  uses  the  edge 
detection  operator,  double  Hires  are  shown  at  lower  iterations,  The  weaker  line  of  the 
double  lines  disappeared  gradually  in  the  lelaxation  process  with  4  successive  iterations 
as  shown  in  Figure  14.  Figure  17  shows  the  results  for  3db  and  6db  of  S/N.  A  desirable 
example  of  edge  detection  is  shown  in  Figure  18. 

C.  MODIFIED  EDGE  DETECTION  ALGORITHM  RXEG 

The  results  of  RXEG  shows  that  this  algorithm  works  as  an  edge  detection  algo¬ 
rithm.  It  can  detect  tlie  edges  in  the  lofargram,  even  though  the  S/N  is  very  low.  Double 
lines  in  Figure  15  (b)  shows  an  example  of  the  edge  detection. 
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Two  iterations(l,2)  of  reinforcement  process  in  RXLN 


Artificial  lofargram  of  6dh  and  its  initial  magnitude  from  RXLN 


Figure  13.  Initial  magnitude  of  Fig  8  (a)  in  RXEG.  The  Prewitt  operator  is  used 
for  initial  probability. 
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(b) 

Figure  14.  Four  iterations  (1,2, 3,4)  of  RXEG 
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(a)  (b) 


Figure  15.  Artificial  lofargram  Mith  6db  S/N  and  its  initial  magnitude 


(c)  (d) 


Figure  16.  Results  of  four  iterations(l.2,3,4)  of  reinforcement  from  Fig  15 
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Figure  17.  Results  of  tuo  tofargrams,  (a)  3db  (b)  6db. 


Results  of  the  reinforceiiient  with  5  iterations.  However,  the  purpose  of  lofargram 
processing  is  to  detect  the  line  of  the  signal  for  the  interpretation  of  the  lofargram.  The 
desirable  ideal  example  of  the  line  detection  is  ohown  in  Figure  19.  To  accomplish  this, 
some  changes  in  RXHG  are  required.  Many  algorithms  were  tested  to  generate  the  edge 
magnitude  and  line  direction  for  the  initial  probability.  The  new  magnitude  and  direction 
obtained  from  EGPR  were  adopted  in  the  modified  algorithm. 


(a)  (b) 

Figure  18.  Example  of  edge  detection  process  (a)  input  (b)  output 


1.  Edge  magnitude  output 

The  gray  level  of  the  edge  enhanced  magnitude  of  each  pixel  are  processed  by 
masks.  Nine  masks  are  formed  around  a  pixel  in  the  shape  of  pentagon,  hexagon  and 
rectangular  as  shown  in  Figure  20.  The  average  values  and  the  variances  within  the 
masks  around  a  pixel  are  calculated.  'I  hen.  the  average  value  of  the  mask  with  minimum 
variance  is  assigned  as  the  gray  level  of  the  center  pixel. 


(a)  (b) 

Figure  19.  Example  of  line  detection  process  (a)  input  (b)  output 
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Figure  20.  Nine  masks  from  EG  PR 


2.  Line  direction  output 

The  initial  line  direction  of  each  pixel  is  determined  from  the  direction  of  the 
mask  with  the  minimum  variance.  Eight  line  directions  can  be  associ  itcd  with  each  im¬ 
age  point.  The  directions  are  symmetric  to  the  origin.  Therefore,  only  four  of  them  are 
considered  as  0’  =  180’ ,  45’  =  225’,  90’  =  270’  and  135’  =  315’.  The  direction  for  an  edge 
is  the  direction  along  the  edge  with  the  dark  side  of  the  edge  on  the  left  and  the  bright 
side  on  the  right  side.  This  means  that  the  edge  direction  0’  and  180’  arc  different  di- 
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rections.  However,  the  directions  for  line  detection  from  masks  in  EGPR  do  not  have 
dark  side  and  bright  side.  Tlie  direction  in  EGPR  simply  implies  the  direction  of  a  line 
without  dark  side  and  bright  side.  This  means  tiiat  the  line  direction  0’  and  180"  are  the 
same  direction  as  shown  in  Figure  21.  The  edge  direction  with  bright  on  the  right  side 
arc  shown  in  Figure  22. 


Figure  21.  An  example  of  line  direction 

A 

dark  /  bright 

/ 

bright  f  dark 

V 

Figure  22.  An  example  of  edge  direction 

I'he  four  directions  generated  from  the  EGPR  are  shown  in  Figure  13.  But,  start  with 
the  second  iteration  after  the  first  rela.xation  process,  infinite  line  directions  are  consid¬ 
ered  at  each  pi.xel. 

3.  Pixel  and  neiglibor  interaction 

Edge  detection  in  RXEG  detects  double  edges  of  the  signal  line  in  the  lofargram. 
It  can  detect  the  edges  reliably  when  the  S;N  is  high.  But,  the  purpose  here  is  to  detect 
thin  line  as  shown  in  Figure  19.  For  the  thin  line  detection,  it  is  necessary  to  construct 
a  compatibility  function  with  proper  compatibility  values  for  the  interactions  between 
different  classes.  Basically,  the  same  compatibility  function  as  that  of  the  RXEG  can  be 
used  for  the  interactions  betw'een  lines.  The  only  difierence  between  the  line  direction 
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90 


0 

Figure  23.  Four  possible  line  direclious 

and  the  edge  direction  is  that  the  directionality  of  a  line  greater  than  180’  is  considered 
as  the  same  as  its  symmetric  value  througlt  tite  origin.  In  other  words,  the  compatibilities 
of  lines  between  tiie  center  pixel  and  the  ncigljboring  pixels  can  be  handled  in  the  same 
way  as  in  RXIiG. 

'I  he  center  pixel  k  could  have  lour  kinds  of  interactions  with  its  neighbor  pixel 
/,  i.c,  linc/linc,  linc/no-iinc,  no-linc/linc  and  no-linc/no-linc. 
rt.  Linclline  interaction 

Let  k  be  the  centered  pixel  position,  and  /  the  neighbor  position.  Let  a  be 
the  line  direction  at  k,  (i  the  line  direction  at  the  neighbor  pixel  /,  y  be  the  direction  of 
line  joining  k  io  I,  D  the  chessboard  distance  from  k  to  /,  i.c,  ntax(  |  k-l^  | ,  1  k,-ly  | )  . 
Then  the  li  le/line  reinforcement  process  between  k  and  /  has  strength  given  by 

c{k,l]l,l)  =  cos(o:-y)  cos(/?-y)/2^  (4. 1 ) 

which  is  similar  to  eq  (3.4).  Here  /  represents  the  line  class  and  n  represents  non-line 
class,  1 0  sec  the  significance  of  this  definition,  a  few  simple  examples  are  considered  in 
Table  4.  Results  changed  with  dilferent  compatibility  value  of  parallel  lines.  In  the  ex¬ 
periment,  when  a  higher  value  of  compatibility  is  chosen  for  parallel  lines,  background 
noise  increases,  while  a  lower  value  is  chosen  when  the  opposite  situation  takes  place. 
b.  Lincjno-line  interaction 

'I  he  line  probability  at  pixel  k  is  weakened  by  the  no-linc  probability  at  pixel 
/  to  the  degree  c{k,l’,l,n)  defined  by 

c{k,l]l,n)  =  minCO,  -  cos(2a-2y)/2^]  (4.2) 

which  is  similar  to  eq  (3.5).  A  few'  examples  are  shown  in  Table  5. 


35 


c.  No-UnelHm  interaction 

I'he  no-iinc  probability  at  pixel  k  is  alTectcd  by  the  line  probability  at  the 
neighbor  pixel  /  to  the  degree  c{k,n\l,t)  denned  by 


c{k,iv,l,l)  = 


(l--cos(2/?-2v)) 


>n+i 


(4.3) 


which  is  similar  to  eq  (3.6).  A  few  examples  are  shown  in  Table  6. 


case 

a 

P 

Y 

cos{a  -  y)  cos(P  -y) 

1 

1 

90 

90 

90 

1/2'^ 

■Hi 

90 

0 

0 

0 

n 

90 

90 

0 

— 

Tabic  4.  Examples  of  line  /  line  interaction 


case 

a 

Y 

-  cos(2a  -  2y)  /  2*^ 

c{k,  1;  1,  n) 

n 

90 

0 

1/2” 

0 

n 

90 

90 

- 1/2'^ 

-1/2'^ 

n 

90 

270 

-1/2”  , 

- 1  /2’^ 

Table  5.  Examples  of  line  /  no-line  interaction 


case 

a 

P 

Y 

(1  -  COS(2P  -  2y))  /  2'^ 

•1 

0 

90 

0 

1/2'^ 

t 

1 

0 

90 

90 

0 

Table  6.  Examples  of  no-Iine  /  line  interaction 
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d.  No-lmcino-linc  interaction 

The  no-Hne  probabilities  at  the  centered  pixel  k  and  that  of  the  neighbor 
pixel  /  reinforce  each  other  to  the  degree  c{k,n\l,n)  defined  by 

c{k,n\l,n)  =  -\  (4.4) 


which  is  similar  to  eq  (3.7) 

4.  Combined  leinfui  cement  process 

For  each  pixel  k,  the  net  effect  of  its  neighboring  pixei..  an  its  probability  p{k,[) 
and  no-linc  probability  p(/<,«)  =  \-p{k,l)  is  computed  as  follows 


I 

I 


(4.5) 


I 


+^Qp'{l,n)c{k,n;l,n)  (4.6) 

I 


where  C,,  C„  Q,  Q  are  constants  who.sc  sum  is  taken  to  be  one.  The  standard  values 
needed  here  are  C,  =  0.866,  Cj  =  0.124,  C,  =  0.005  ,  and  Q  =  0.005.  'fhe  results  of  the 
iteration  process  are  somewhat  sensitive  to  the  choice  of  the  C's.  For  examjile,  if  C,  is 
too  large,  the  line  will  thicken  and  will  be  extended  into  no-linc  points;  while  if  Q  is  too 
large,  gaps  will  appear  at  weak  spots  in  the  lines  and  at  sharp  angles., 


<iW)  = 


I  I  + 1  I 


(4.7) 


q'{k,n) 


I  7V<.0  I  +  I  ?'(^r'0  I 


(4.8) 
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p'(k,0=p(kMiHW)] 

(4.9) 

p'(k,»)  =  p'(k,n)[l+q'(k,n);] 

(4.10) 

r+I-.  ,  /(^.O 

^  p\kMp\k,n) 

(4.11) 

This  process  is  then  iterated  with  replacing  and  l-p'*'{l<,f)  replacing 

p'{li,n).  We  also  compute  the  estimated  line  direction  at  each  point 

cos(0(/()) 

+^p\lMkAl,f)cosm)  (4.12) 

i 

A'y{k)  -  U'p\k,l)  sm{0{k)) 

+^/(  WA, /;/,/)  sin(0(0)  (4.13) 

I 

0^+'(/<)  =  tan-'(A',(/<)/A',(/<))  (4.14) 

where  iV  is  a  constant.  For  large  values  of  the  constant  W,  is  close  to  0;  while  for 
small  values,  it  is  strongly  influenced  by  the  neighboring  0's. 

5.  Experimental  results 

Figure  24  shows  the  initial  magnitude  of  the  input  image  of  Figure  8  (a).  Figure 
25  shows  the  results  for  3db  and  6db  of  S/N.  I'igure  26  shows  the  results  of  four  iter¬ 
ations  of  the  reinforcement  applied  to  Figure  24  using  the  modilied  RXEG.  Figure  27 
shows  the  artificial  lofargram  with  6db  of  S/N  and  its  initial  magnitude.  Figure  28  shows 
the  results  of  four  iterations  of  reinforcement  with  6db  S/N  artificial  lofargram.  Com¬ 
paring  Figure  13  and  Figure  24  reveals  that  the  initial  magnitudes  of  the  modified  RXliG 
are  improved  by  using  the  EGI’R  line  detection.  The  initial  magnitude  of  modified 
RXFG  shows  stronger  signal  line  in  Figure  24  when  compared  to  Figure  13.  However, 
the  end  results  of  reinforcement  do  not  improve  significantly.  From  Figure  28,  it  is  ob¬ 
vious  that  a  double  line  still  exists  and  that  the  lines  are  thicker  than  those  of  the  RXEG 


algorithm.  In  this  experiment,  the  purpose  is  to  obtain  single  thin  line  from  the 
lofargram.  But  the  result  is  not  fully  satisfactory.  A  few  things  could  be  dis  issed  with 
respect  to  the  suspected  reasons  of  the  double  lines  and  of  the  thicker  line  in  the  result. 
The  initial  magnitude  as  shown  in  Figure  27  (b)  was  caused  by  the  Prewitt  operator  of 
the  RXEG.  In  Figure  1.5  the  result  of  RXFG  thows  strongly  the  splitting  line.  Therefore, 
the  reinforcement  process  couldn't  be  the  problem  of  causing  double  lines.  Figure  29 
shows  the  diagram  of  the  modified  RXEG.  As  for  the  possible  cause  of  the  thicker  line, 
the  modified  initial  magnitude  may  be  the  reasot..  Comparison  of  Fig  15  (b)  and  Fig  27 
(b)  reveals  the  improvement  of  the  initial  magnitude.  The  initial  magnitude  of  the  mod¬ 
ified  RXEG  shows  stronger  and  thicker  lines.  Consequently,  the  results  of  the  re¬ 
inforcement  at  a  later  stage  applied  to  this  initial  magnitude  shov;*'  thicker  lines. 
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Figure  24.  Inilial  magnitude  of  Fig  15  (a)  In  modified  RXEG 


To  solve  the  double  line  problem,  it  is  necessary  to  find  a  method  to  yield  the 
initial  magnitude  in  single  lines.  Two  dilfeicnt  method  were  tested  for  this  purpose.  In 
the  first  method  each  pixel  was  divided  by  the  maximum  value  over  the  entire  picture 
as  follow. 
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(4.15) 


rikJn) 

max  p{k,ln) 

k 

This  method  removed  noise  very  quickly.  The  initial  magnitudes  at  each  pixel  are  cal¬ 
culated  by  dividing  each  pixel  by  the  maximum  values  over  the  entire  picture.  Therefore, 
very  small  values  are  assigned  to  noise  and  the  high  values  are  assigned  to  the  signal  line 
when  the  S/N  was  high.  On  the  other  hand,  when  S/N  is  low,  the  high  values  are  as¬ 
signed  to  both  the  noise  and  the  signal  because  the  difference  of  gray  level  between  the 
noise  and  the  signal  is  small.  As  a  result,  this  method  is  still  not  suitable  to  lofargram 
processing.  In  the  second  method,  each  pixel  was  divided  by  the  maximum  value  of  the 
local  neighborhood  as  follow 


pjkjn) 

maxp{ljn) 


(4.16) 


Initial  magnitude  in  this  method  depends  on  the  local  maximum  value.  In  the  exper¬ 
iments  it  is  realized  the  initial  magnitude  could  be  white  or  dark  entirely  depending  on 
the  characteristic  of  the  noise  of  the  lofargram.  This  second  method  was  not  appropriate 
for  the  calculation  of  the  initial  magnitude  before  the  relaxation  process. 

In  this  chapter  various  results  of  the  algorithm  studied  were  discussed.  Thick 
line  problem  were  discovered.  Several  attempts  has  been  made  to  solve  this  problem. 
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(d) 

Fioure  26.  Four  iterations  (1,2, 3, 4)  of  modified  RXEG  (continue) 
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J 


(a)  (b) 

Figure  27.  Arlificial  lofargram  >^ith  6cib  S/N  and  its  initial  magnitude 


(c)  (d) 

Figure  28.  Results  of  four  iterations  (1,2, 3, 4)  from  Fig  27  in  modified  RXEG 


Figure  29.  The  procedure  of  niudified  RXEG 
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V.  CONCLUSIONS  AND  RECOMMENDATION 


The  fact  that  each  algorithm  has  its  own  peculiar  characteristics  in  the  lofargrani 
processing  was  found  in  this  study.  1  hese  characteristics  arc  very  much  dependent  on  the 
S/N  of  the  input  lofargrani.  Algorithms  RXLN  and  RXLAl  using  similar  ideas  can  de¬ 
tect  lines  in  the  loHirgram  when  the  S/N  is  relatively  high.  But,  the  result  were  very  dis¬ 
couraging  when  the  S/N  was  low  and  noise  smoothing  was  done  slowly  in  the  iterations 
of  the  reinforcement.  Our  main  concern  of  the  lofargrani  processing  is  to  detect  a  thin 
line  when  S/N  is  low.  ‘Ihcrcforc,  thc.se  algorithms  arc  not  suitable  for  lofargrani  process. 
On  the  other  hand,  the  RXilG  algorithm  results  in  good  lofargram  processing.  The  only 
shott  coming  is  that  it  is  an  edge  detection  algorithm.  It  smoothes  the  noise  and  en¬ 
hances  the  line  components  in  the  lofaigram  even  at  low  S/N.  The  last  algorithm 
studied  in  this  thesis  is  the  inodincd  RXIid  algorithm.  Several  method  were  tested  to 
solve  the  double  line  problem.  But,  double  line  problem  still  e,\ists  in  the  modified  algo¬ 
rithm.  ‘1 0  solve  this  problem,  it  is  necessary  to  search  for  an  algorithm  that  can  yield  thin 
single  line  in  the  image  of  the  initial  magnitude. 
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APPENDIX .  THE  PROGRAM  OF  MODIFIED  RXEG 

INTEX3ERM  ISX/256/,  ISY/256/,  ITER/4/,  ISNX/5/,ISNy/5/,  lSW/1/ 
INTBGERM  IP  ( 256 , 256 ) ,  IPRB(  256, 256 ) ,  JP  ( 256 , 256 ) 

DIMENSION  C(4),CTYPE{4,2),THET0(256,256> 

DIMENSION  PRBEl (256,256), PRBE2( 256 , 256 ) , 0( 8 ) 

DIMENSION  THETl ( 256 , 256 ) ,THET2( 256 , 256 ) 

BYTE  G(256) 

REAL  LENGTH,  H/3 . 0/,MAX,MIN 

OPEN ( UNIT-1, NAME- ' BEAMOO . BIN ' , TYPE- ' OLD ' , ACCESS- ' DIRECT ' , 

2  RECORDSIZE-64,MAXREC-256) 

OPEN(UNIT-2,NAME-*TLOFAR20-4X5.DAT‘  ,TYPE-'NEW'  ,ACCESS-' DIRECT’  , 
2  RECORDSIZE-64,MAXREC-256) 


DO  10  1-1,256 
READ(1'I)G 
DO  20  J-1,256 

IP(I,J)-G(J) 

1F(IP(I,J).LT.0)  IP(I,J>-IP{l,J)+256 
20  CONTINUE 

10  CONTINUE 

CALL  RXEG  ( IP ,  THETl ,  THET2 ,  PRBEl ,  PRBE2 , 3 SX ,  ISY ,  ISNX ,  ISNY ,  1  TER , 
&  ISM, TEMP, JP,THET0) 

MAX-PRBE2(1,1) 

MIN-PRBE2(1,1) 

DO  30  1-1,256 
DO  40  J-1,256 

IF(MAX.LT.PRBE2(I,J))  MAX-PRBE2( 1 , J> 

IF (MIN , GT . PRBE2 ( 1 , J ) )  MIN-PRBE2 ( 1 , J ) 

40  CONTINUE 

30  CONTINUE 

LENGTH- (MAX-MIN) 

DO  50  1-1,256 
DO  60  J-1,256 

IPRB ( I , J ) -JNINT ( ( ( PRSF2 ( I ,  J  ) -MIN )/LENGTH ) * 24 3 . ) 

60  CONTINUE 

50  CONTINUE 

DO  70  1-1,256 
DO  80  J-1,256 

IF(IPRB(I,J),GT.127)  THEN 
G(J)-IPRB(I,J)-25e 
ELSE 

G(J)-IPRB(I,J) 

ENDIF 

80  CONTINUE 

WRITE(2'I)G 
70  CONTINUE 

END 

SUBROUTINE  RXEG (IP , THETl , THET2 , PRBEl , PRBE2 , ISX , ISY 
S  , ISNX, ISNY, ITER, ISW, TEMP, JP,THET0) 

C 

C 

CP  Edge  reinforcement  by  releocation  method. 
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c 

CS  CALL  RXBG(IP,THBT1,THET2,PRBE1,PRBE2,ISX,ISY 

CS  /1SNX,ISNY, ITER, ISW, TEMP) 

C 

CK  RELAXATION,  EDGE 


IP(ISX,ISY)  ! 

Input  image 

(I) 

THBT1(ISX,ISY)! 

Gradient  direction  after  ITER-1 
iterations 

(W) 

THi:,T2(ISX,ISY)j 

Gradient  direction  after  ITER  iterations 

(0) 

PRBE1(ISX,ISY): 

Edge  probability  for  each  pixel  after 
lTER-1  iterations 

(W) 

PRBE2(ISX,ISY)8 

Edge  pr^>ability  for  each  pixel  after 
ITER-1  iterations 

(0) 

ISNX,ISNY  ; 

Size  of  neighborhood  being  considered 
(odd  number,  for  example:  5,5) 

(I) 

ITER  : 

Number  of  iterations 

(I) 

ISH  t 

Switch  (1  -  standard,  2  -  noise  reduction 
type 

(I) 

Reference 

(31!  B.  J.  Schachter,  A.  Lev,  S  W.  Zuc)cer,  and  A.  Rosenfeld,  "An 
application  of  relaxation  methods  to  edge  reinforcement," 

IEEE  Trans.  Syst.,  Man,  Cybern.,  vol.  SMC-7,  pp.  813-816, 

Nov.  1977. 

JULY  1979.  PROGRAMMED  BY  K.SAKAUE 

(1)  When  computing  the  Initial  values  for  edge  probabifty  and 
edge  direction  using  an  operator  other  than  the  Prewitt 
operator,  rewrite  lower- level  subroutine  RXEP.  Any  operators 
that  ca  obtain  a  differential  (Dx,Dy)  can  be  used. 

(2)  The  processing  characteristic  largely  depends  upon  selection 
of  coefficients  Cl,  through  C4  (See  the  reference.)  The  two 
standard  )clngs  of  values  have  already  been  set  in  array  CTYPE 
us.'.ng  data  statements,  and  values  can  be  changed  with  ISW. 
vr.'  n  using  other  values  or  changing  values  depending  upon 
ti.t  number  of  iteration  tiroes,  rewrite  the  program  so  a  proper 
value  is  set  to  the  argxunent  C  in  lower- level  subroutine  RXEI. 
The  following  condition  however,  must  be  met. 

C(l)+C(2)+C(3)+C(4)-1 

(3)  The  coefficient  W  used  when  updating  (THETl)  is  set  in  data 
statement  to  3.0.  The  same  value  was  also  used  in  the  ref  feme 
This  value  may  be  changed  as  neccessary. 

DIMENSION  IP(ISX,ISY),JP(ISX,ISY) 

DIMENSION  THETl ( ISX , ISY ) , THETO ( ISX , ISY ) 

DIMENSION  THET2(ISX,ISY) 

DIMENSION  PRBE1(ISX,ISY) 

DIMENSION  PRBE2 ( ISX , ISY ) , TEMP ( 256 , 256 ) 

DIMENSION  CTYPE (4, 2) 

DIMENSION  C(4) 

DATA  CTYPE/. 866, .124, .005,. 005 

4  ,.706,. 176,. 059,. 059/ 

DATA  W/3.0/ 

CF  INITIAL  EDGE  VALUES  AND  INITIAL  PROBABILITIES 
C 
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CALL  BGPR(IP,JP,1SX,ISY,THET1,THET0) 

CALL  RXro(lP,THETl,PRBEl,lSX,lSY,TEMP) 

C 

CP  RELAXATION  UPDATING  PROCESS 
C 

AMN-0. 

AMX-0. 

DO  9  1-1,4 

C(I)-CTYPE(I,ISW) 

9  CONTINUE 
DO  10  I-1,ITER 

IF(I.BQ.l)GO  TO  11 
DO  12  IY-1,1SY 
DO  12  1X-1,1SX 

THETl < IX , lY ) -THET2 ( IX , lY) 

PRBEl ( IX , lY ) -PRBE2 ( IX , lY ) 

12  CONTINUE 
11  CONTINUE 

CALL  RXEI  ( THETl ,  THET2 ,  PRBEl ,  PRBE2 ,  ISX ,  ISY ,  ISNX ,  ISNY , C ,  W) 
10  CONTINUE 
RETURN 
END 


C 

C 

CS 

C 

CP 

C 

CK 

C 

CA 

CA 


SUBROUTINE  EGPR ( IP , JP , ISX, ISY , THETl , THETO ) 
Copyright  (c)  1983  by  AIST  MITI(JAPAN) 

CALL  EGPR(IP,JP,1SX,ISY) 

Edge  preserving  smoothing  operation. 

EDGE  PRESERVING  SHOOTHING,  SMOOTHING 

IP(ISX,ISY)  :  Input  image  array 

JP(ISX,ISY)  :  Output  image  array 


(I> 

(0) 


M 

CM 

CM 

C 


C 


The  program  is  considerably  long,  but  the  main  part  consists 
of  about  120  lines  in  the  first  half  and  the  remainder  is  used 
as  routines  for  edge  processing. 

DIMENSION  IP( ISX, ISY) , THETl (ISX, ISY) , THETO ( ISX,  ISY) 

DIMENSION  JP(I3X,ISY) 

DIMENSION  A(9),  V(9),0(9) 

DIMENSION  K(5,5),  L(5,5) 


EQUIVALENCE 

& 

& 

EQUIVALENCE 

& 

& 


(A(1),A1),  (A(2),A2),  (A(3),A3), 
(A(4),A4),  (A(5),A5),  (A(6),A6), 
(A(7),A7),  (A(8),A8),  (A(9),A9) 
(V(1),V1),  (V(2),V2),  (V(3),V3), 
(V(4),V4),  (V(5),V5>,  (V(6),V6), 
<V(7},V7),  (V{8),V8),  (V(9),V9) 


C 


EQUIVALENCE 

&  {K11,K(1,1)>,{K21,K(2,1)),(K31,K(3,1)),(K41,K(4,1)),(K51,K(5,1)), 


t  (K14,K(1,4)> 
t  (K15,K{1,5)) 
EQUIVALENCE 
&  (L11,L(1,1)) 
4  (L12,L(1,2)) 
&  (L13,L(1,3)) 
6  (L14,L(1,4)) 
&  {L15,L(1,5)) 


/{K24,K(2,4) 

.(K25,K(2,5) 


),<K34,K<3,4)> 

),(K35,K(3,5)) 


,<K44,K(4,4> 

,(K45,K(4,5) 


>,(K54,K(5,4>) 

),(K55,K(5,5)> 


,(L21,L<2,1) 

,(L22,L(2,2) 

,(L23,L(2,3) 

,(L24,L(2,4) 

,(L25,L(2,5) 


),(L31,L(3,1)) 

),<L32,L(3,2)) 

),(L33,L(3,3)) 

),{L34,L<3,4)) 

),(L35,L<3,5)) 


,(L41,L(4,1) 

,(L42/L(4,2) 

,(L43,L(4,3) 

,(L44/L(4,4) 

,(L45/L(4,5) 


),(L51,L(5,1)> 

),(L52,L(5,2)) 

),(L53,L(5,3)) 

),(L54,L(5,4)) 

),(L55,L(5,5)) 


C 

CF 


AVE1(K1,K2,K3,K4,K5,K6,K7,K8,K9)  - 
(  FLOAT(Kl+K2+K3+K4+K5+K6+K7+K8+K9 )/9 . 0 

AVE2(K1,K2,K3,K4,K5,K6,K7)  - 
«  FLOAT(K1+K2+K3+K4+K5+K6+K7)/7.0 

VAR1(A1,L1,L2,L3,L4,L5,L6,L7,L8,L9)  - 
«  7 . 0* (FLOAT(Ll+L2+L3+L4+L5+L6+L7+Le+L9 )-Al*Al*9 . 0 )/9 . 0 

VAR2(A1,L1,L2,L3,L4,L5,L6,L7)  - 
&  FLOAT(L1+L2+L3+L4+L5+L6+L7)-A1*A1*7.0 


*****  Initialization 


ISXO 

ISYO 

ISXl 

ISYl 

ISX2 

ISY2 

ISX3 

ISY3 


ISX 
ISY 
ISXO-1 
ISYO-1 
ISXO-2 
ISYO- 2 
ISXO-3 
ISYO- 3 


C 

CF 


*****  Smoothing  process  for  the  region  (3,3>- 
DO  60  IY-3,ISY2 
IYM2  -  IY-2 
IYP2  -  IY+2 
DO  60  IX-3,ISX2 
IXM2  -  IX-2 
IXP2  -  IX+2 


■(ISX-2,ISY-2) 


KY  -  0 

DO  20  JY-IYM2,IYP2 
KY  -  KY+1 
KX  -  0 

DO  20  JX-IXM2,IXP2 
KX  -  KX+1 
IPD  -  IP(JX/JY) 
K(KX,KY)  -  IPD 
L(KX,KY)  -  IPD*IPD 
20  CONTINUE 


C 

CF 


***  Calculation  of  mean  values  and  variances  for  nine  masks 


A1  -  AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2  -  AVE2(K21,K31,K41,K22,K32,K42,K33) 

A3  -  AVE2(K41,K51,K32,K42,K52,K33,K43) 

A4  -  AVE2(K42,K52,K33/K43,K53,K44,K54) 

A5  -  AVE2(K33,K43,K34/K44,K54,K45,K55) 

A6  -  AVE2(K33,K24,K34,K44,K25,K35,K45) 

A7  -  AVE2(K23,K33,K14,K24,K34,K15,K25) 

A8  -  AVE2(K12,K22,K13,K23,K33,K14,K24) 

A9  -  AVE2(Kll,K21,K12,K22,K32,K23,K33) 

VI  -  VAR1(A1,L22,L32,L42,L23,L33,L43,L24;L34,L44) 


V2  -  VAR2(A2,L21,L31,M1,L22,L32,L42,L33) 

V3  -  VAR2(A3,L41,L51,L32,L42,L52,L33,L43) 

V4  -  VAR2(A4,L.42,L52,L33,L43,L53,I/44,L54) 

V5  -  VAR2(A5,L33,L43,L34,L44,L54,L45/L55) 

V6  -  VAR2(A6,L33,L24,L34,L44,L25,L35,M5) 

V7  -  VAR2(A7,L23,L33,U4,L24,L34/L15,L25) 

V8  -  VAR2(A8,L12,L22,L13,L23,L33,L14,L24) 

V9  -  VAR2(A9,L11,L21,L12,L22,L32,L23,L33) 

CF  ***  SET  THE  DIRECTION  FOR  NINE  MASKS 
Q(2)-3.14 
Q(3)-3.925 
Q(4)-4.71 
Q(5)-5.495 
Q(6)-0.0 
Q(7)-0.785 
Q(8)-1.57 
Q(9)-2.355 

CF  ***  Mean  gray  value  of  the  mask  with  minimum  variance  — >  JP(IX/IY) 
RMIN  -  VI 
MI  -  1 
DO  40  1-2,9 

IF  (RMIN  .LE.  V(I))  GO  TO  40 
RMIN  -  V(I) 

MI  -  I 

40  CONTINUE 

IF(MI.E0.6)THEN 
THET0(IX,IY)-0.0 
GO  TO  45 

ELSE  IF(MI.EQ.7)THEN 
THET0(IX,IY)-0.785 
GO  TO  45 

ELSE  IF(MI.E0.8)THEN 
THET0(IX,IY)-1.57 
GO  TO  45 

ELSE  IF(MI.EQ.9)THEN 
THET0(IX,IY)-2.355 
GO  TO  45 

ELSE  IF(MI.E0.2)THEN 
THET0(IX,IY)-3.14 
GO  TO  45 

ELSE  IF (MI. EQ. 3) THEN 
THET0(IX,lY)-3.925 
GO  TO  45 

ELSE  IF (MI. EQ. 4) THEN 
THET0(IX,IY)-4.71 
GO  TO  45 

ELSE  IF (MI. EQ. 5) THEN 
THET0(IX,IY)-5.495 
GO  TO  45 

ELSE  IF(MI.EQ.1>THEN 
GO  TO  41 
END  IF 

41  CONTINUE 
SMALL-V(2) 

ASMALL-Q(2) 

DO  42  1-3,9 
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IF(SMALL.GT.V(1) )THEN 
SMALJ>V(  I ) 


ASMALL-O(I) 

ENDIF 
42  CONTINUE 

THETO(IX,IY)-ASMALL 
45  CONTINUE 

JP(IX,IY)  -  A(MI)+0.5 


60  CONTINUE 
C 

CF  *****  Operation  for  the  four  corners 
C 

C - 

C 

DO  80  IY-1;4 
DO  80  IX-1,4 
IPD  -  IP(IX,1Y) 

K(IX,IY)  -  IPD 
L(IX,IY)  -  IPD* IPD 
80  CONTINUE 
C 

CF*****  Jp(l,l)  ***** 

JP(1,1)  -  AVE2(Kll,K21,K12,K22,K32,K23,K33)+0.5 
C 

CF*****  JP(2,1)  ***** 

A1  -  AVE2(K21,K31,K22,K32,K42,K33,K43) 

A2  -  AVE2(K21,K12,K22,K32,K13,K23,K33) 

VI  -  VAR2(A1,L21,L31,L22,L32,L42,L33,L43) 

V2  -  VAR2(A2,L21,L12,L22,L32,L13,L23,L33) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(2,1)  -  RMIN+0.5 
C 

CF*****  JP(1,2)  ***** 

A1  -  AVE2(K21,K31,K12,K22,K32,K23,K33) 

A2  -  AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI  -  VAR2(A1,L21,L31,L12,L22,L32,L23,L33) 

V2  -  VAR2(A2,L12,L22,L13,L23,L33,L24,L34) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(1,2)  -  RMIN+0,5 
C 

CF*****  JP(2,2)  ***** 

A1  -  AVE1(K11,K21,K31,K12,K22,K32,K13,K23,K33) 

A2  -  AVE2(K31,K41,K22,K32,K42,K33,K43> 

A3  -  AVE2(K22,K32,K23,K33,K43,K34,K44) 

A4  -  AVE2(K22,K13,K23,K33,K14,K24,K34) 

VI  -  VAR1(A1,L11,L21,L31,L12,L22,L32,L13,L23,L33) 
V2  -  VAR2(A2,L31,L41,L22,L32,L42,L33,L43) 

V3  -  VAR2(A3,L22,L32,L23,L33,L43,L34,L44) 

V4  -  VAR2(A4,L22,L13,L23,L33,L14,L24,L34) 

C 

RMIN  -  VI 
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non 


MI  -  1 

DO  100  1-2,4 

IF  (RMIN  .LE.  V(I))  GO  TO  100 
RMIN  -  V<I) 

MI  -  I 

100  CONTINUE 

JP{2,2)  -  A(MI)+0.5 


DO  120  n-1/4 
KX  -  0 

DO  120  IX-ISX3,ISX0 

KX  -  KX+1 

IPD  -  1P(IX,IY) 

K(KX,IY)  -  IPD 
L(KX,IY)  -  IPD* IPD 
120  CONTINUE 
C 

CF*****  JP(ISX-1,1)  ***** 

A1  -  AVE2<K31,K22,K32,K42,K23,K33#K43) 

A2  -  AVE2(K21,K31,K12,K22,K32,Kp»X23) 

VI  -  VAR2(A1,L31,L22,L32,M2,L23,L33,W3) 

V2  -  VAR2(A2,Ii21,L31,L12,L22,L32,L13,L23) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP{ISXl,l)  -  RMIN+0.5 

c 

rv*****  JP(ISX”1,2)  ***** 

A1  -  AVE1(K21,K31,K41,K22,R32,K42,K23,K33,K43) 

A2  -  AVE2(K32,K23,K33,K43,K24,K34,K44) 

A3  -  AVE2(K22,K32,K13,K23,K33,K14,K24) 

A4  -  AVE2(K11,K21,K12,K22,K32/K13,K23) 

VI  -  VAR1(A1,L21,L31,L41,L22,L32,L42,L23,L33,L43) 
V2  -  VAR2(A2,L32,L23,L33,L43,L24,L34,M4) 

V3  -  VAR2(A3,L22,L32,L13,L23,L33,L14,L24) 

V4  -  VAR2{A4 ,L11,L21,L12,L22,L32,L13,L23) 

C 

RMIN  -  VI 
Ml  -  I 

DO  140  1-2,4 

IF  (RMIN  .LE.  V(I))  GO  TO  140 
RMIN  -  V(I) 

MI  -  1 

140  CONTINUE 
C 

JP(ISX1,2)  -  A(Ml)+0.5 

-’’aTC2;K32,M2!k23.1!33,K43,K24,M4) 

A2  -  AVE2(K21,K31,K22,K32,K42,K23,Kp) 

VI  -  VAR2(A1,L32,L42,L23,L33,L43,L24,L34) 

V2  -  VAR2(A2,L21,L31,L22,L32,L42,L23,L33) 


uuo 


RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(ISX0,2>  -  RMIN+0.5 


KY  -  0 

DO  160  IY-ISY3,ISY0 
KY  -  KY+1 
DO  160  IX-1,4 
IPD  -  1P(IX,1Y) 

K(IX,KY)  -  IPD 
L(IX,KY)  -  IPDMPD 
160  CONTINUE 
C 

CF*****  JP(1,ISY-1)  ***** 

A1  -  AVE2(K21,K31,K12,K22,K32,K13,K23) 

A2  -  AVE2(K22,K32,K13,K23,K33,K24,K34) . 

VI  -  VAR2(A1,L21/L31,L12,L22,L32,L13,L23) 

V2  -  VAR2(A2,L22,L32,L13,L23,L33/L24,L34) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(1,ISY1)  -  RMIN+0.5 
C 

CF*****  JP(2,ISY-1)  ***** 

hi  -  AVE1(K12,K22,K32,K13,K23,K33,K14,K24,K34) 

A2  -  AVE2(K11,K21,K31,K12,K22,K32,K23) 

A3  -  AVE2(K31,K41,K22,K32,K42,K23,K33) 

A4  -  AVE2(K32,K42,K23,K33,K43,K34,K44) 

VI  -  VAR1(A1,L12,L22,L32,L13,L23,L33/L14,L24,L34). 
V2  -  VAR2(A2,L11,L21,L31,L12,L22,L32,L23> 

V3  -  VAR2(A3,L31,L41,L22,L32,L42,L23,L33) 

V4  -  VAR2(A4,L32,L42,L23,L33,L43,L34,L44) 

C 

RMIN  -  VI 
MI  -  1 
DO  180  1-2,4 

IF  (RMIN  .LE.  V(I))  GO  TO  180 
RMIN  -  V(I) 

MI  -  I 

180  CONTINUE 
C 

JP(2,ISY1)  -  A(MI)+0.5 
C 

CF*****  JP(1,ISY)  ***** 

JP(1,ISY)  -  AVE2(K22,K32,K13,K23,K33,K14,K24)+0.5 
C 

CF*****  ■JP(2,ISY)  ***** 

A1  -  AVE2(K12,K22,K32,K13,K23,K33,K24) 

A2  -  AVE2(K32,K42,K23,K33,K43,K24,K34) 

VI  -  VAR2(A1,L12,L22,L32,L13,L23,L33,L24) 

V2  -  VAR2(A2,L32,L42,L23,L33,L43,L24,L34) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(2,ISY0)  -  RMIN+0.5 
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KY  -  0 

DO  200  IY-1SY3,ISY0 
KY  -  KY+1 
KX  -  0 

DO  200  1X-ISX3,ISX0 
'  KX  -  KX+1 

IPD  -  IP(IX,IY) 

K(KX,KY)  -  IPD 
L(KX,KY)  -  IPD*IPD 
200  CONTINUE 
C 

CF*****  JP(ISX-1,1SY-1)  ***** 

A1  -  AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2  -  AVE2(K21,K31,K41,K22,K32,K42,K33) 

A3  -  AVE2(K12,K22,K13,K23,K33,K14,K24) 

A4  -  AVE2(K11,K21,K12,K22,K32,K23,K33) 

VI  -  VAR1(A1,L22,L32,L42,L23,L33,L43,L24,L34,L44) 

V2  -  VAR2(A2,L21,L31,L41,L22,L32,L42,L33) 

V3  -  VAR2(A3,L12,L22,L13,L23,L33,L14,L24) 

V4  -  VAR2(A4,L11,L21,L12,L22,L32,L23,L33) 

C 

RMIN  -  VI 
MI  -  I 
DO  220  1-2,4 

IF  (RMIN  .LE.  V(l))  GO  TO  220 
RMIN  -  V(I) 

MI  -  1 

220  CONTINUE 
C 

JP(ISX1,ISY1)  -  A(MI)+0.5 
C 

CF*****  JP(ISX,ISY-1)  ***** 

A1  -  AVE2(K22,K32,K23,K33,K43,K24,K34) 

A2  -  AVE2(K21,K31,K22,K32,K42,K33,K43) 

VI  -  VAR2(A1,L22,L32,L23,L33,L43,L24,L34) 

V2  -  VAR2(A2,L21,L31,L22,L32,L42,L33,L43) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(ISX0,ISY1)  -  RMIN+0.5 
C 

CF*****  JP(ISX-1,ISY)  ***** 

A1  -  AVE2(K22,K32,K42,K23,K33,K43,K34) 

A2  -  AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI  -  VAR2(A1,L22,L32,L42,L23,L33,L43,L34) 

V2  -  VAR2(A2,L12,L22,L13,L23,L33,L24,L34) 

C 

RMIN  -  A1 

IF  (VI  .GT.  V2)  RMIN  -  A2 
JP(ISX1,ISY0)  -  RMIN+0.5 
C 

CF*****  JP(1SX,ISY)  ***** 

JP(ISX0,ISY0)  -  AVE2(K22,K32,K23,K33,K43,K34,K44)+0.5 
C 

CF  *****  Marginal  operation 
C 

C - 

C 
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DO  360  1X-3,ISX2 
IXM2  -  IX-2 
IXP2  -  IX+2 
C 

Q - 

C 

DO  240  JY-1,4 
KX  -  0 

DO  240  JX-IXM2,IXP2 

KX  -  KX+1 

IPD  -  1P(JX,JY) 

K(KX,JY)  -  IPD 
L(KX,JY)  -  IPD*IPD 
240  CONTINUE 
C 

CF*****  JP(1X,1)  ***** 

A1  -  AVE2(K31,K41,K32,K42,K52,K43,K53) 

A2  -  AVE2(K31,K22,K32,K42,K23,K33,K43) 

A3  -  AVE2(K21,K31,K12,K22,K32,K13,K23)  . 

VI  -  VAR2(A1,L31,I/41,L32,L42,L52/L43,L53) 

V2  -  VAR2<A2,L31,L22,L32,L42,L23,L33,L43) 

V3  -  VAR2(A3,L21,L31,L12,L22,L32,L13,L23) 

C 

RMIN  -  VI 
MI  -  1 
DO  260  1-2,3 

IF  (RMIN  .LE.  V(I))  GO  TO  260 
RMIN  -  V(I) 

MI  -  I 

260  CONTINUE 
C 

JP(IX,1)  -  A(MI)+0.5 
C 

CF*****  JP(IX,2)  ***** 

A1  -  AVE1(K21,K31,K41,K22,K32,K42,K23,K33,K43) 

A2  -  AVE2(K41,K51,K32,K42,K52,K43,K53) 

A3  -  AVE2(K32,K42,K33,K43,K53,K44,K54) 

A4  -  AVE2(K32,K23,K33,K43,K24,K34,K44) 

A5  -  AVE2(K22,K32,K13,K23,K33,K14,F24> 

A6  -  AVE2(K11,K21,K12,K22,K32,K13,K23) 

VI  -  VAR1(A1,L21,L31,L41,L22,L32,L42,L23,L33,L43) 
V2  -  VAR2(A2,L41,L51,L32,L42,L52,L43,L53) 

V3  -  VAR2(A3,L32,L42,L33,L43,L53,L44,L54) 

V4  -  VAR2(A4,L32,L23,L33,L43,L24,L34,L44) 

V5  -  VAR2(A5,L22,L32,L13,L23,L33,L14,L24) 

V6  -  VAR2(A6,L11,L21,L12,L22,L32,L13,L23) 

C 

RMIN  -  VI 
MI  -  1 
DO  280  1-2,6 

IF  (RMIN  .LE.  V(I))  GO  TO  280 
RMIN  -  V(I) 

MI  -  I 

280  CONTINUE 
C 

JP(IX,2)  -  A(MI)+0.5 
C 

C - 
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C 

KY  -  0 

DO  300  JY-ISY3,ISY0 
KY  -  KY+1 
KX  -  0 

DO  300  JX-IXM2,IXP2  ‘ 

KX  -  KX+1 
IPD  -  IP(JX,JY) 

K(KX,KY)  -  IPD 
L(KX,KY)  -  IPD*IPD 
300  CONTINUE 
C 

CT*****  JP(IX,ISY-1)  ***** 

M  -  AVE1(K22,K32,K42,K23,K33,K43,K24,K34,K44) 

A2  -  AVE2(K42,KS2,K33«K43,K53,K44,KS4) 

A3  -  AVE2(K41,K51,K32,K42,K52,K33,K43) 

A4  -  AVE2(K21,K31,K41,K22,K32,K42,K33) 

A5  -  AVE2(K11,K21,K12/K22,K32,K23,K33) 

A6  -  AVE2(K12,K22,K13/K23,K33,K14,K24) 

VI  -  VARl(Al,L22,li32/L42,L23,L33,l/43,L24,L34,L44> 
V2  -  VAR2(A2»L42,L52,L33,L43,L53,L44,L54) 

V3  -  VAR2(A3,L41,L51,L32,L42,L52,L33,L43) 

V4  -  VAR2(A4,L21,L31,L41,L22,L32,L42,L33) 

V5  -  VAR2(A5,L11,L21,L12,L22,L32,L23,L33) 

V6  -  VAR2(A6,L12,L22,L13,L23,L33,L14,L24) 

C 

RMIN  -  VI 
MI  -  1 
DO  320  1-2,6 

IF  (RMIN  .LE.  V(I))  GO  TO  320 
RMIN  -  V(I) 

MI  -  I 

320  CONTINUE 
C 

JP(IX,ISY1)  -  A(MI)+0.5 
C 

CF*****  JP(IX,ISY)  ***** 

A1  -  AVE2(K42,K52,K33,K43,K53,K34,K44) 

A2  -  AVE2(K22,K32,K42,K23,K33,K43,K34) 

A3  -  AVE2(K12,K22,K13,K23,K33,K24,K34) 

VI  -  VAR2(A1,L42,L52,L33,L43,L53,L34,L44) 

V2  -  VAR2(A2,L22,L32,L42,L23,L33,L43,L34) 

V3  -  VAR2{A3,L12,L22,L13,L23,L33,L24,L34> 

C 

RMIN  -  VI 
MI  -  1 
DO  340  1-2,3 

IF  (RMIN  .LE.  V(I))  GO  TO  340 
RMIN  -  V(I) 

MI  -  I 

340  CONTINUE 
C 

JP(IX,ISY0)  -  A(MI)+0.5 
C 

360  CONTINUE 
C 
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DO  500  1Y-3,ISY2 
IYM2  -  lY-2 
IYP2  -  IY+2 
C 

Q - 

C 

KY  -  0 

DO  380  JY-IYM2,IYP2 
KY  -  KY+1 
DO  380  JX-1,4 
IPD  -  IP(JX,JY) 

K(JX,Ky)  -  IPD 
L(JX,KY)  -  IPD*IPD 
380  CONTINUE 
C 

CF****«  JP(1,IY)  ***** 

A1  -  AVE2(K21,K31,K12,K22,K32,K13,K23) 

A2  -  AVE2(K22,K32,K13,K23,K33,K24,K34> 

A3  -  AVE2(K13,K23,K14,K24,K34,K25,K35) 

VI  -  VAR2(A1,L21,L31,L12,L22,L32,L13,L23) 

V2  -  VAR2(A2,L22,L32,L13,L23/L33,L24,L34) 

V3  -  VAR2(A3,L13,L23,L14,L24,L34,L25,L35) 

C 

RMIN  -  VI 
MI  -  1 
DO  400  1-2,3 

IF  (RMIN  .Lt.  V(I))  GO  TO  400 
RMIN  -  V(I) 

MI  -  I 

400  CONTINUE 
C 

JP(1,IY)  -  A(MI)+0.5 
C 

CF*****  JP(2,IY)  ***** 

A1  -  AVE1(K12,K22,K32,K13,K23,K33,K14,K24,K34) 

A2  -  AVE2(K11,K21,K31,K12,K22,K32,K23) 

A3  -  AVE2(K31,K41,K22,K32,K42,K23,K33) 

A4  -  AVE2(K32,K42,K23,K33,K43,K34,K44) 

A5  -  AVE2(K23,K33,K24,K34,K44,K35,M5) 

A6  -  AVE2(K23,K14,K24,K34,K15,K25,K35) 

VI  -  VAR1(A1,L12,L22,L32,L13,L23,L33,L14,L24,L34) 
V2  -  VAR2(A2,L11,L21,L31,L12,L22,L32,L23) 

V3  -  VAR2(A3,L31,L41,L22,L32,L42,L23,L33) 

V4  -  VAR2(A4,L32,L42,L23,L33,L43,L34,L44) 

V5  -  VAR2(A5,L23,L33,L24,L34,L44,L35,L45) 

V6  -  VAR2(A6,L23,L14,L24,L34,L15,L25,L35) 

C 

RMIN  -  VI 
MI  -  1 
DO  420  1-2,6 

IF  (RMIN  .LE.  V(I))  GO  TO  420 
RMIN  -  V(I) 

MI  -  I 

420  CONTINUE 
C 

JP(2,iy)  -  A(MI)+0.5 
C 

C - 
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c 

KY  -  0 

DO  440  JY-IYM2,IYP2 
KY  -  KY+1 
KX  -  0 

DO  440  JX-ISX3,ISX0 

KX  -  KX+1 

IPD  -  IP{JX,JY) 

K(KX,KY)  -  IPD 
L(KX,KY)  -  IPD*1PD 
440  CONTINUE 
C 

CF*****  JP(ISX-1,IY)  ***** 

A1  -  AVE1(K22,K32/K42,K23,K33,K43,K24,K34,K44^ 

A2  -  AVE2(K21«K31,K41,R22,K33,K42,K33) 

A3  -  AVE2(K11,K21,K12,K22,K32,K23,K33) 

A4  -  AVE2(K12,K22/K13,K23,X33,K14,X24) 

A5  -  AVE2(K23,K33,K14,K24,K34,K15,K25) 

A6  -  AVE2(K33,K24,K34,K44,K25,K35,K45) 

VI  -  VAR1{A1,L22,L32,L42,L23,L33,L43,L24,L34,L44) 
V2  -  VAR2(A2,L21,L31,L41,L22,L32,L42,L33) 

V3  -  VAR2(A3,L11,L21,L12,L22,L32,L23,L33) 

V4  -  VAR2(A4,L12,L22,L13,L23,L33,L14,L24) 

V5  -  VAR2(A5,L23,L33,L14,L24,L34,L15,L25) 

V6  -  VAR2(A6,L33,L24,L34,L44,L25,L35,L45) 

C 

RMIN  -  VI 
MI  -  1 
DO  460  1-2,6 

IF  (RMIN  .LE.  V(I))  GO  TO  460 
RMIN  -  V(I) 

MI  -  I 

460  CONTINUE 
C 

JP(ISX1,IY)  -  A(MI)+0.5 
C 

CF*****  JP(ISX,1Y)  ***** 

A1  -  AVE2(K21,K31,K22,K32,K42,K33,K43) 

A2  -  AVE2(K22,K32,K23,K33,K43,K24,K34) 

A3  -  AVE2(K33,K43,K24,K34,K44,K25,K35) 

VI  -  VAR2(Al,L2l,L31,L22,L32,L42,L33,L43) 

V2  -  VAR2(A2,L22,L32,L23,L33,L43,L24,L34) 

V3  -  VAR2(A3,L33,L43,L24,L34,L44,L25,L35) 

C 

RMIN  -  VI 
MI  -  1 
DO  480  1-2,3 

IF  (RMIN  .LE.  V(I))  GO  TO  480 
RMIN  -  V(I) 

MI  -  I 

480  CONTINUE 
C 

JP(ISX0,IY)  -  A(MI)+0.5 
C 

500  CONTINUE 

DO  600  1-1,256 
DO  600  J-1,256 
IP(I,J)-JP(I,J) 

600  CONTINUE 
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DO  601  IX-1,2 
DO  601  IY-1,256 
THET0(IX,IY)-0.0 

601  CONTINUE 

DO  602  IX-1,256 
DO  602  IY-1,2 

THETO(IX,IY)-0.0 

602  CONTINUE 

DO  603  IX-255,256 
DO  603  IY-1,256 
THET0(IX,IY)-0.0 

603  CONTINUE 

DO  604  IX-1,256 
DO  604  IY-255,256 
THETO(IX,IY)-0.0 

604  CONTINUE 

DO  605  IX-1,256 
DO  605  IY-1,256 

THETl ( IX , lY ) -THETO  ( IX , lY ) 

605  CONTINUE 

DO  606  IX-1,256 
DO  606  IY-1,256 

IF{THET1( IX, lY)  .GE.  3 . 14 )THEN 

THETl ( IX , lY ) -THETl ( IX , I Y ) - 3 . 14 
ELSE 

THETl ( IX , lY ) -THETl ( IX, lY  > 

ENDIF 

606  CONTINUE 
C 

CF  *****  RETURN 
RETURN 
END 

SUBROUTINE  RXEP( IP, THETl, PRBEl, ISX, ISY,TEMP) 

C 

C  Copyright  (c)  1983  by  AIST  MITI (JAPAN) 

C 

CP  Computes  initial  parameter  and  initial  probability  for  edge 
CP  reinforcement  by  relaxation  method  (Lower-level  routine 

CP  for  subroutine  RXEG). 

C 

CS  CALL  RXEP(IP, THETl, PRBEl, ISX, ISY) 

C 

CK  RELAXATION,  EDGE 

C 

CA  IP(ISX,ISY)  !  Input  image 

CA  THETl ( ISX, ISY) :  Initial  value  of  gradient  direction  for 

CA  each  pixel 

CA  PRBE1(ISX,ISY) s  Initial  value  of  edge  probability  for 

CA  each  pixel 

C 

CN  Reference 

CN  @1!  B.  J.  Schachter,  A.  Lev,  S.  W.  ZucJcer,  and  A.  Rosen f eld,  "An 
CN  application  of  relaxation  methods  to  edge  reinforcement," 

CN  IEEE  Trans.  Syst.,  Man,  Cybern. ,  vol.  SMC-7,  pp.  813-816, 

CN  Nov.  1977. 

C 

CD  JULY  1979.  PROGRAMMED  BY  K.SAKAUE 

C 

CM  A  small  real  number  -1E30  is  set  in  variable  AMX  as  an  initial 
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CM  value.  When  using  a  computer  with  different  word  lengths/  care 

CM  should  be  taken. 

C 

DIMENSION  1P(ISX,ISY)/TEMP(ISX,ISY) 

DIMENSION  THET1(ISX,ISY) 

DIMENSION  PRBE1(ISX/ISY) 

DIMENSION  NE(3/3) 

C 

CF  INITIAL  EDGE  VALUES 

CF  BY  PREWITT  OPERATOR  (DIFFERENTIAL  TYPE) 

C 

AMX— 1E30 
DO  10  IY-1,ISY 
DO  10  IX-1,ISX 
DO  11  J-1,3 
DO  11  1-1,3 
Jl-IY+J-2 
ll-IX+I-2 

IF(J1.LT.1.0R.J1.GT.1SY)G0  TO  12 
IF(I1.LT.1.0R.I1.GT.ISX)G0  TO  12 
NE(I,J)-1P(11/J1) 

11  CONTINUE 

DX-FLOAT(NE(3/l)+NE(3,2)+NE(3,3)-NE(l,l)-NE(l/2)-NE(l/3)) 
DY-FLOAT(NE( 1 / 3 )+NE( 2 , 3 )+NE( 3 , 3 )-NE(  1 , 1 )-NE( 2 / 1)-NE( 3 , 1) ) 
AMG-SQRT ( DX*  *  2+DY  *  *  2 ) 

I F ( AMG . GT . AMX ) AMX- AMG 
PRBE1(1X/IY)-AMG 

IF(DX.LE.0E0. AND. DY.LE.OEO)DX-0. 0001 


GO  TO  10 
12  CONTINUE 

PRBE1(IX/IY)-0.001 


10 


20 


C 

C 

C 

CS 

C 

CP 

CP 

CP 

C 

CK 

C 

CA 

CA 

CA 

CA 


CONTINUE 
DO  20  IY-1,ISY 
DO  20  IX-1,ISX 

PRBEl ( IX , IY)-PRBE1 ( IX, IY)/AMX 
PRBE1(IX,IY)-AMAX1( .01,AMIN1( . 9,PRBE1( IX, lY) ) ) 

TTMP ( IX , lY ) -PRBEl ( IX , lY ) 

CONTINUE 

RETURN 

END 

SUBROUTINE  RXEI  ( THETl ,  THET2 ,  PRBEl ,  PRBE2 ,  ISX ,  ISY ,  ISNX ,  ISNY ,  C ,  W ) 


Copyright  (c)  1983  by  AIST  MITI( JAPAN) 


CALL  RXEI ( THETl , THET2 , PRBEl , PRBE2 , ISX , ISY , ISNX , ISNY , C ,W) 

Updates  edge  parameter  and  edge  probed>ility  for  edge 
reinforcement  by  relaxation  method. 

(Lower-level  routine  for  subroutine  RXEG) 


RELAXATION,  EDGE 


THET1(ISX,ISY) 

THET2(ISX,ISY) 

PRBE1(ISX,ISY) 

PRBE2(ISX,ISY) 


Gradient  direction 
Gradient  direction 
Edge  probabilities 
Edge  probabilities 


(before  updating) 
(after  updating) 
(before  updating) 
(after  updating) 
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C 


c 

c 

c 


c 

CF 

C 


C 

CF 

CF 

c 


ISNX/ISNY  !  Size  of  neighborhood  being  considered 

(limited  to  odd  numl'er  three  or  higher) 
C(4)  :  Parameter  Cl  through  C4  (See  the 

reference. ) 

.866, .124, .005, .005  (standard) 

.706,  .176,  .059  (noise  cleaning) 

W  :  Parameter  W  (See  the  reference.) 

3.0  for  standard  use. 

DIMENSION  THET1(ISX,ISY) 

DIMENSION  THET2(ISX,-TSY) 

DIMENSION  PRBE1(1SX,ISY) 

DIMENSION  PRBE2(ISX,ISY) 

DIMENSION  C(4) 

DATA  HPAI/1. 570796327/ 


ISFTX-ISNX/2+1 
ISFTY-ISNY/2+1 
DO  10  IY-1,ISY 
DO  10  IX-1,ISX 

ALPHA-THET1(IX,IY) 

PXY-PRBE1(IX, lY) 

Ql-0.0 

Q2-0.0 

DHX-W*  PXY*C0S ( ALPHA ) 

DHY-W*PXY*SIN( ALPHA) 

DO  11  J-1,ISNY 
DO  11  I-1,ISNX 
lU-IX+l-ISFTX 
IV-IY+J-ISFTY 

IF(rV.EQ.lY.AND.lU.EO.IX)GO  TO  11 

IF(IV.LT.1.0R.IV.GT.ISY)G0  TO  11 

IF(IU.LT.1.0R.IU.GT.ISX)G0  TO  11 

IDY-IV-IY 

IDX-IU-IX 

JI>-IABS(IDY) 

ID-IABS(IDX) 

LD-MAX0(ID,JD) 

LD  !  CHESSBOARD  DISTANCE  FROM  (IX,IY)  TO  (IU,R') 

ALD-FLOAT(LD) 

D2-1.0/2.0**ALD 

BETA-THET1(IU,IV) 

DY-FLOAT(IDY) 

DX-FLOAT( IDX) 

IF(IDX.EQ.0. AND. IDY.EQ.0)DX-. 00001 
GAMMA«ATAN2(DY,DX) 

PUV-PRBE1(IU,IV) 

COMPATIBILITY  COEFFICIENTS 
EDGE/EDGE  INTERACTION 

IF ( ALPHA . EQ . BETA . AND . BETA . EQ . GAMMA ) THEN 
REE- 1.0 


ELSE  IF ( ALPHA .  EQ .  BETA .  AND .  GAMMA ,  NE .  BETA )  THEN 
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REE-0 . 5 
ELSE 

REE-COS  ( ALPHA-GAMMA)  *COS  ( BETA-GAMMA )  *  D2 
ENDir 

c 

CF  INTERACTION  WITH  NONEDGES 

C 

AG2-2 . 0* (ALPHA-GAMMA) 

AG2— COS(AG2)*D2 
REN-AMIN1(0.0,AG2) 

BG2-2 . 0* (BETA-GAMMA) 

RNE- ( 1 . O-COS ( BG2 ) ) *D2*0 . 5 
RNN-D2 
C 

CF  COMBINED  REINFORCEMENT  PROCESS 

C 

Q1-Q1+C( 1 ) *PUV*REE+C( 2) * ( 1 . 0-PUV) *REN 
02-Q2+C ( 3 ) * PUV*RNE+C( 4 ) * ( 1 . 0-PUV) *RNN 
PREE-PUV*REE 
DHX-DHX+PREE*COS ( BETA ) 

DHY-DHY+PREE*  SIN ( BETA) 

11  CONTINUE 

00-ABS(Q1)+ABS(02) 

Ql-Ql/QQ 

Q2-Q2/QO 

P1D-PXY*(1.0+Q1) 

P2D-(1.0-PXY)*(1.0+Q2) 

PNEW-PlD/( P1D+P2D) 

IF(DHX.LE.0E0. AND. DHY.LE.OEO)DHY-0. 00001 
THNEW-ATAN2 ( DHY , DHX ) 

PRBE2(IX,IY)-PNEW 
THET2(IX,1Y)-THNEW 
10  CONTINUE 
RETURN 
END 
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